Karlsruher Institut für Technologie

**Schriftenreihe Kontinuumsmechanik im Maschinenbau**

# Jannick Kuhn

**25**

J. Kuhn **Micromechanical modeling of polycrystalline metals**

Microstructure modeling and crystal plasticity parameter identification for predicting the cyclic mechanical behavior of polycrystalline metals

25

Jannick Kuhn

#### **Microstructure modeling and crystal plasticity parameter identification for predicting the cyclic mechanical behavior of polycrystalline metals**

#### **Schriftenreihe Kontinuumsmechanik im Maschinenbau Band 25**

Karlsruher Institut für Technologie (KIT) Institut für Technische Mechanik Bereich Kontinuumsmechanik

Hrsg. Prof. Dr.-Ing. habil. Thomas Böhlke

Eine Übersicht aller bisher in dieser Schriftenreihe erschienenen Bände finden Sie am Ende des Buchs.

## **Microstructure modeling and crystal plasticity parameter identification for predicting the cyclic mechanical behavior of polycrystalline metals**

by Jannick Kuhn

Karlsruher Institut für Technologie Institut für Technische Mechanik Bereich Kontinuumsmechanik

Microstructure modeling and crystal plasticity parameter identification for predicting the cyclic mechanical behavior of polycrystalline metals

Zur Erlangung des akademischen Grades eines Doktor-Ingenieurs von der KIT-Fakultät für Maschinenbau des Karlsruher Instituts für Technologie (KIT) genehmigte Dissertation

von Jannick Kuhn, M.Sc.

Tag der mündlichen Prüfung: 3. November 2022 Hauptreferent: Jun.-Prof. Dr. rer. nat. Matti Schneider Korreferent: Prof. Dr.-Ing. Sebastian Münstermann

**Impressum**

Karlsruher Institut für Technologie (KIT) KIT Scientific Publishing Straße am Forum 2 D-76131 Karlsruhe

KIT Scientific Publishing is a registered trademark of Karlsruhe Institute of Technology. Reprint using the book cover is not allowed.

www.ksp.kit.edu

*This document – excluding parts marked otherwise, the cover, pictures and graphs – is licensed under a Creative Commons Attribution-Share Alike 4.0 International License (CC BY-SA 4.0): https://creativecommons.org/licenses/by-sa/4.0/deed.en*

*The cover page is licensed under a Creative Commons Attribution-No Derivatives 4.0 International License (CC BY-ND 4.0): https://creativecommons.org/licenses/by-nd/4.0/deed.en*

Print on Demand 2023 – Gedruckt auf FSC-zertifiziertem Papier

ISSN 2192-693X ISBN 978-3-7315-1272-1 DOI 10.5445/KSP/1000154640

## **Zusammenfassung**

Aufgrund ihrer Kombination aus Festigkeit und Zähigkeit werden polykristalline Metalle oft für Strukturbauteile verwendet, welche anspruchsvollen Betriebsbedingungen, z. B. zyklischen mechanischen Belastungen, ausgesetzt sind. Wegen der hohen Sicherheitsvorgaben sowie der ökologischen und wirtschaftlichen Anforderungen ist es unerlässlich für eine sichere und kosteneffiziente Konstruktion, die eingesetzten Werkstoffe für die jeweiligen Belastungsbedingungen zu optimieren. Insbesondere bei zyklischer mechanischer Beanspruchung werden die Eigenschaften des Werkstoffes stark durch die zugrunde liegende Mikrostruktur beeinflusst. Darum muss diese Wechselwirkung bei der Gestaltung einer Komponente berücksichtigt werden. Die numerische Homogenisierung bietet eine Möglichkeit den Einfluss der Mikrostruktur auf die makroskopischen mechanischen Eigenschaften von polykristallinen Werkstoffen abzubilden. Die vorliegende Arbeit beschäftigt sich mit der Erzeugung von Abbildungen polykristalliner Mikrostrukturen welche als Eingangsgröße für die numerische Homogenisierung dienen können. Des Weiteren wird für die robuste und effiziente Parametrisierung des zugrunde liegenden mikromechanischen Modells anhand makroskopischer Experimente ein Optimierungsansatz vorgeschlagen, welcher auf Methoden des maschinellen Lernens beruht.

Die Abbildung der Mikrostruktur, auch synthetische Mikrostruktur genannt, muss relevante Eigenschaften, wie zum Beispiel die Korngrößenverteilung und kristallografische Textur des realen Polykristalls widerspiegeln. Hierfür werden in dieser Arbeit sogenannte Laguerre-Tessellierungen verwendet, welche die Domäne der Mikrostruktur in

einzelnen Zellen zerlegt. Die Unterteilung erfolgt auf Basis von Nuklei, welche zuvor in der Domäne zufällig verteilt wurden und den Tessellierungsgewichten. Durch das Lösen eines konvexen Optimierungsproblems können die Gewichte so bestimmt werden, dass die Zellen eine vorgeschriebene Größe haben. Hierdurch können synthetische Mikrostrukturen erzeugt werden, welche experimentell bestimmte Korngrößenverteilungen abbilden. Ebenso können mehrphasige synthetische Polykristalle mit vorgegebenen Volumenanteilen erzeugt werden. Moderne, gradientenbasierte Löser werden im Hinblick auf die Anzahl der Iteration und die Gesamtlaufzeit, welche sie zur Lösung des Optimierungsproblems benötigen, miteinander verglichen.

Im nächsten Schritt wird eine Methodik entwickelt um den berechneten Laguerre-Zellen kristallografische Orientierungen zuzuweisen. Hierzu wird eine Kostenfunktion definiert, welche den Unterschied zwischen vorgegebenen und aus der synthetischen Mikrostruktur berechneten Koeffizienten einer Fourier-Reihenentwicklung der kristallographischen Orientierungsverteilungsfunktion beschreibt. Zur Minimierung des Unterschiedes zwischen diesen sogenannten Texturkoeffizienten wird ein Gradientenabstiegsverfahren verwendet, welches iterativ die kristallographischen Orientierungen der synthetischen Mikrostrukturen anpasst. Mit diesem Ansatz wird die Möglichkeit untersucht, isotropes sowie anisotropes mechanisches Verhalten zu reproduzieren, sowohl bei zugrundeliegenden gleich- als auch log-normal-verteilten Korngrößen. Aufgrund der optimierten Orientierungen reproduziert die synthetische Mikrostruktur das gewünschte mechanische Verhalten mit vergleichsweise wenigen Körnern, d.h. kleinen Volumenelementen.

Um auf Basis der zuvor generierten synthetischen Mikrostrukturen das mechanische Verhalten vorherzusagen, muss das verwendete mikromechanische Materialmodell parametrisiert werden. Es wird vorgeschlagen mithilfe von makroskopischen polykristallinen Experimenten die mikromechanischen Parameter zu identifizieren. Hieraus ergibt sich ein inverses Optimierungsproblem, wobei die Auswertung der sich daraus ergebende Zielfunktion ressourcenintensiv ist. Zur Minimierung der nötigen Funktionsaufrufe wird die Bayes'sche Optimierungsstrategie vorgeschlagen, welche ein statistisch informiertes Ersatzmodell der Kostenfunktion aufbaut und mithilfe geeigneter Erfassungsfunktionen die Erforschung und Ausnutzung ausbalanciert. Die Leistungsfähigkeit und Robustheit der Bayes'schen Strategie wird untersucht und mit anderen, für ähnliche Probleme verwendete, Optimierungsalgorithmen verglichen.

## **Summary**

Due to their combination of strength and toughness, polycrystalline metals are often used for structural components that are subjected to demanding operating conditions, e.g. cyclic mechanical loads. Owing to the high safety requirements as well as environmental and economic restrictions, tailoring these materials to their respective loading conditions for a safe and cost-efficient design seems imperative. Especially in the case of cyclic mechanical loading, the behavior of the material is strongly influenced by the underlying microstructure. Thus, this interaction must be taken into account when designing a component. Computational homogenization offers a way to capture the influence of the microstructure on the macroscopic mechanical properties of polycrystalline materials. In this work we address generating representations of polycrystalline microstructures which can be used as input for computational homogenization. Furthermore, we propose using an optimization approach based on machine-learning, for a robust and efficient parameterization of the underlying micromechanical model, based on macroscopic experiments.

The representation of the microstructure, also called synthetic microstructure, must reflect relevant properties, e.g., grain size distribution and crystallographic texture, of the real polycrystal. For this purpose we use so-called Laguerre tessellations, which decomposes the synthetic microstructure domain into individual cells. This decomposition uses, previously in the domain randomly distributed, nuclei and tessellation weights. By solving a convex optimization problem, the weights can be determined so that the cells have a predefined size. This allows generating synthetic microstructures, that reproduce experimentally observed grain size distributions. Similarly, multi-phase synthetic polycrystals with given volume fractions can be generated. We compare different modern gradient-based solvers in terms of the number of iterations and overall run time they need to solve the optimization problem.

In a next step, we develop a methodology to assign crystallographic orientations to the Laguerre cells. For this purpose, we define a cost function, describing the difference between given and from the synthetic microstructure computed coefficients of a Fourier series expansion of the crystallographic orientation distribution function. For minimizing the difference between these, so-called texture coefficients, we use a gradient descent method, which iteratively adjusts the crystallographic orientations of the synthetic microstructures. We investigate the capability of the method to reproduce isotropic as well as anistropic mechanical behavior, both with underlying uniform and log-normal grain size distributions. Due to the optimized orientations, the synthetic microstructure reproduces the desired mechanical behavior with few grains, i.e., small volume elements.

To predict the mechanical behavior based on the previously generated synthetic microstructures, the micromechanical material model must be parameterized. We propose using macroscopic polycrystalline experiments to identify the micromechanical parameters. This results in an inverse optimization problem with an objective function, which is resource intensive to evaluate. We propose to use a Bayesian optimization strategy in order to minimize the necessary function evaluations, which builds a statistically informed surrogate model of the cost function and uses appropriate acquisition functions to balance exploration and exploitation. We investigate the performance and robustness of the Bayesian optimization approach and compare it to other optimization algorithms used for similar problems.

# **Acknowledgments**

My sincerest thanks belongs to Matti Schneider for his continuous guidance and support. I am very grateful for all the discussions and conversations, whether on current research topics or personal matters. I am grateful to Sebastian Münstermann for his interest in this work, for taking over the role as a second examiner and his complementary view on micromechanics. I want to thank Christoph Kirchlechner for chairing the thesis committee and the material scientific perspective he provided. Last but not least, I would like to express my gratitude towards Thomas Böhlke for the insights into the mechanics of polycrystalline metals and providing me the opportunity to pursue my thesis at the Institute of Engineering Mechanics (ITM).

A big thank you goes to Petra Sonnweber-Ribic for her unwavering support during my time at the Robert Bosch GmbH. Her confidence in my work, the encouragement and the support in all organizational matters were extremely important to me.

I would like to thank all of the colleagues at the Center of Research and Advance Engineering of the Robert Bosch GmbH for numerous fruitful discussions and a stimulating working atmosphere. Many thanks to Lars Epple for all the discussions about Electron Backscatter Diffraction methods, Elke Bachmeyer for supporting me with the metallographic analysis and Thomas Waldenmaier for conducting numerous heat treatments. I want to say thank you to Fabian Welschinger for the discussions about FeelMath and the HPC cluster.

My thanks to all the colleagues of AMP3 for deep insights into modern engineering metals.

For the many on-site but also virtual discussions and the great after-work exchanges I would like to thank my fellow doctoral candidates at the Robert Bosch GmbH, Benjamin Schäfer, Nikolai Arnaudov, Aleksander Eric, Marius Graf, Marius Wolf, Argha Protim Dey and Erik Natkowski.

Towards the colleagues at the ITM I would like to say thank you for lively coffee breaks and the fun activities outside of work. I really enjoyed my time as an external PhD candidate at the ITM thanks to Hannes Erdle, Alexander Dyck, Sebastian Gajek, Johannes Gisy, Johannes Görthofer, Patrick Hessmann, Tobias Karl, Johannes Kreusten, Maximilian Kraus, Celine Lauff, Jonas Lendvai, Alok Mehta, Lennart Risthaus, Daniel Wicht, Andreas Prahs and Loredana Kehrer. A special thanks to Felix Ernesti for being an extraordinary colleague and an even better friend. For all the help with organizational matters I would like to thank Tom-Alexander Langhoff, Helga Betsarkis, Ariane van Elst and Ute Schlumberger-Maas.

A giant thank you to all my friends, who supported me during this time and helped me to not get lost in my studies.

I am eternally grateful towards my parents for always encouraging me to pursue my dreams and who supported me throughout my life.

With all my heart I would like to thank my wife Hannah for her support, patience and love. I could always count on you to have my back.

Karlsruhe, December 2022 Jannick Kuhn

# **Contents**




# **Chapter 1 Introduction**

## **1.1 Motivation**

Meeting the high safety, environmental and customer demands in the automotive industry can be quite challenging. Especially, since the loading conditions can be complex and environmental restrictions require a light-weight design in order to reduce energy consumption. To construct components fulfilling these restrictions, suitable materials have to be identified.

Low-alloyed, high-strength carbon steels are often employed because of their high strength and toughness while being cost-effective, see Horvath (2021) for an overview on applications in the automotive sector. Designing components made from these materials in a *safe* way, while taking into account the need for weight reduction, requires considering the loading conditions which the component has to face. Approximately 90% of metallic components fail due to fatigue (Richard et al., 2014), which describes failure under cyclic loading conditions *below* the static failure limit. To prevent fatigue upon everyday use, it is necessary to assess the cyclic mechanical behavior, including the fatigue limit, of the component under consideration.

Typically, the fatigue limit is determined via mechanical experiments. A pertinent engineering tool to assess the fatigue properties of a component are Wöhler diagrams, which show the applied strain or stress

amplitude over the endured cycles. For the same design and material of a component, there is a stochastic influence on lifetimes, i.e., for the same loading conditions, the cycles to failure scatter in the Wöhler diagram. The cause for this scattering lies in the inhomogeneous microstructure of a component's material, which dictates the deformation behavior under given loading conditions, see Sec. 2.2 for a discussion about polycrystalline metals and François (1989) as well as Krupp (2007) for further insights.

Thus, designing metallic components to endure cyclic loading conditions requires taking the underlying microstructure into account. Fig. 1.1 shows an image of a ferritic microstructure (Arnaudov et al., 2020), obtained by Electron Back Scatter Diffraction (EBSD) methods.

**Figure 1.1:** Image of a polycrystalline microstructure obtained by Electron Backscatter Diffraction (EBSD) (Arnaudov et al., 2020) and postprocessed by MTex (Nolze and Hielscher, 2016).

Each crystallite is colored according to its orientation with respect to the laboratory system. We observe, even for a single phase polycrystal, a spatial distribution of crystallographic orientations within the

microstructure. As the mechanical properties depend on the crystallite orientation, there is a spatial distribution of mechanical properties within the polycrystal, i.e., a spatial variation of the deformation behavior. This heterogeneity in the microstructure leads to the scatter typically observed in fatigue tests (Krupp, 2007).

To reduce time- and cost-intensive experiments, required to capture the stochastic influence on fatigue caused by the microstructure, analytical and computational homogenization techniques may be used. These methods allow determining the effective mechanical behavior based on available microstructure data. Whereas being computational more demanding than analytical approaches, numerical full-field homogenization schemes provide access to local fields, e.g., stresses and strains as well as plastic deformations, forming in the microstructure. These fields serve as input to predict the fatigue behavior of a material under given loading (Schäfer et al., 2019a;b; Natkowski et al., 2021a;b).

The above methods require input data in terms of a computational domain, representing the microstructure, as well as a material model, describing the deformation behavior on the microscale. This thesis provides novel tools to supply this input for polycrystalline materials to computational homogenization software, i.e., methods to create computational cells and identify parameters for the micromechanical material model.

## **1.2 Outline**

In Chap. 2 we introduce the foundations of this work. Starting from the basics of continuum mechanics, we introduce a pertinent material model to capture the described deformation behavior. Additionally we provide an overview of the basics of homogenization and the necessary input for computational homogenization. Throughout Chap. 3 to Chap. 5 we introduce the novelties forming the basis of this work.

In Chap. 3 we build upon the work by Bourne et al. (2020). They enabled the fast computation of Laguerre tessellations with cells of prescribed volume fraction by solving a convex optimization problem. We harness the potential of modern gradient-type-solvers to speed up solving this convex problem and compare them for a variety of applications. In addition, we use Anderson acceleration to reduce the computation time when generating so-called centroidal tessellations, i.e., tessellations where the seed of a cell coincides with its centroid. The combination of modern solvers and Anderson acceleration enables the generation of (centroidal) polycrystalline microstructures of industrial complexity within seconds up to a few minutes, improving upon other reported polycrystalline microstructure generators. Indeed, we show that with this method it is possible to generate a polycrystalline morphology which matches experimentally observed statistics, e.g., grain size distribution, with little computational effort.

Extending the previously proposed method to generate tessellations, representing polycrystalline microstructures, we propose a novel method to equip the resulting cells with suitable *crystallographic* orientations in Chap. 4. Using the texture coefficients, i.e., the tensorial coefficients of a Fourier series expansion of the crystallite orientation distribution function, we formulate a cost function which computes the difference between prescribed and realized texture coefficients. To minimize this difference, we use a gradient descent scheme, iteratively adjusting the crystallographic orientations within the synthetic microstructure. We

compare the proposed "*T*exture coefficient *O*ptimization for *P*rescribing orientations" (TOP) method to state-of-the-art approaches. To do so, we compare the effective linear-elastic and the non-linear plastic behavior for a number of realizations and a varying number of grains as well as for the uniform and a textured case. Last but not least we investigate the capability of the proposed method to reproduce a desired mechanical behavior with an additional underlying grain size distribution.

Chap. 5 applies the Bayesian optimization approach, typically used in machine learning, to the problem of identifying the crystal plasticity parameters based on polycrystalline experiments. Because extending the homogenization software by automatic differentiation strategies is difficult (or close to impossible when using commercial software), solving this inverse optimization problem requires minimizing a blackbox function. The Bayesian optimization framework iteratively builds up a statistically informed surrogate model of the underlying function and uses this information to drive the optimization procedure. Testing this approach on different search spaces, we gain insight into the performance of Bayesian optimization and the energy landscape of the underlying black-box function. By comparing the Bayesian approach to commonly used algorithms in materials science as well as performant derivative-free algorithms, we show the efficiency and robustness of the Bayesian optimization when identifying crystal plasticity parameters. The methods and main results are summarized in Chap. 6 with a brief overview of possible applications. As a closing point, future potential fields of research are discussed.

# **Chapter 2 Fundamental concepts**

## **2.1 Basics of continuum mechanics**

### **2.1.1 Objectives**

To predict the cyclic behavior of metallic components we need a suitable tool set, enabling the description of deformation under arbitrary loading conditions. The framework of continuum mechanics provides exactly these tools, allowing the rigorous description of body movement, deformation behavior and their thermodynamic basis. This section is based on Sect. 1 and 2 of (Šilhavý, 1997; Haupt, 2013), Sect. 2 and 3 in (Bertram, 2012) and the work by Truesdell and Noll (2004), without any claim of originality or completeness.

#### **2.1.2 Kinematics**

Kinematics refers to the smooth movement of a material body ℬ ⊂ R in time through a -dimensional Euclidean space (in this thesis we consider = 3). We parameterize the current position of a body ℬ at time by the function : ℬ<sup>0</sup> × [0*,* ] → R 

$$x = \chi(X, t),\tag{2.1}$$

where ∈ R denotes an arbitrary material point position in reference configuration ℬ0. Using the spatial inverse mapping −1

$$X = \chi^{-1}(x, t) \tag{2.2}$$

allows identifying the reference position based on the current placement. The displacement field : ℬ<sup>0</sup> × [0*,* ] → R , defined by

$$u(X,t) = \chi(X,t) - X \tag{2.3}$$

describes the difference between the current and the reference placement, see Fig. 2.1.

**Figure 2.1:** Motion and displacement of a continuum body ℬ.

Describing a field quantity Λ in terms of the current placement, denoted by Λ(*,* ), is called the Eulerian description. Using the reference configuration to describe the same quantity, i.e., Λ(*,* ), is known as the Lagrangian description.

These are connected through (Haupt, 2013)

$$
\Lambda\_E(x, t) = \Lambda\_L(\chi^{-1}(x, t), t), \tag{2.4}
$$

$$
\Lambda\_L(X, t) = \Lambda\_E(\chi(X, t), t). \tag{2.5}
$$

Following Truesdell and Noll (2004), the derivative of Λ with respect to time yields, for the Lagrangian description,

$$
\dot{\Lambda}\_L(X,t) = \frac{\partial \Lambda\_L}{\partial t}(X,t),
\tag{2.6}
$$

and, in the Eulerian description,

$$
\dot{\Lambda}\_E(x,t) = \frac{\partial \Lambda\_E}{\partial t}(x,t) + \frac{\partial \Lambda\_E}{\partial x}(x,t) \cdot \dot{x},\tag{2.7}
$$

where ˙ = (*,* ) denotes the velocity of the material point. For the reader's convenience we will drop the subscripts and in the following and indicate the description by the arguments (*,* ) and (*,* ), respectively.

The deformation gradient : ℬ<sup>0</sup> × [0*,* ] → R × , defined by

$$F(X,t) = \frac{\partial \chi}{\partial X}(X,t),\tag{2.8}$$

is a tensorial field quantity, which describes the variation of deformation in the neighborhood of a material point from reference to current configuration. Thus, it maps the infinitesimal line elements in the reference placement ℬ<sup>0</sup> to the corresponding line elements in the current configuration ℬ

$$dx = F \cdot dX.\tag{2.9}$$

Analogous explicit formulas for the conversion of both area and volume elements undergoing the deformation are available

$$da = \det(F) F^{-T} \cdot dA \tag{2.10}$$

$$dv = \det(F)dV.\tag{2.11}$$

For a rigid body movement in three dimensions without rotation, the corresponding deformation gradient computes to = Id, where Id denotes the identity on rank two tensors.

We quantify the deformation in an infinitesimal neighborhood around a material point , e.g., the change of length and angles between line elements, by the Green's strain (Haupt, 2013)

$$E^{G} = \frac{1}{2} \left( F^{T} F - \text{Id} \right). \tag{2.12}$$

Using the displacement gradient : ℬ<sup>0</sup> × [0*,* ] → R × , defined by

$$H(X,t) = \frac{\partial u}{\partial X}(X,t),\tag{2.13}$$

Green's strain reads

$$E^G = \frac{1}{2} \left( H + H^T + H^T H \right). \tag{2.14}$$

If the displacements are small w.r.t. the typical dimensions in the continuum, i.e.,

$$\mathcal{B}\_t \approx \mathcal{B}\_0 \tag{2.15}$$

the material and spatial coordinates and are approximately equal, i.e., ≈ . Thus, it is not necessary to distinguish between the Lagrangian and Eulerian description, i.e.,

$$
\Lambda(X, t) \approx \Lambda(x, t). \tag{2.16}
$$

For this case, the Frobenius norm ||·|| for all = (*,* ) is small (Bertram, 2012)

$$||H|| = \sqrt{\text{tr}(HH^T)} \ll 1.\tag{2.17}$$

Linearizing the Green's strain around the rest state = 0 gives rise to the infinitesimal strain tensor (Haupt, 2013)

$$
\varepsilon = \frac{1}{2}(H + H^T). \tag{2.18}
$$

#### **2.1.3 Balance equations**

With suitable tools to describe the motion of a continuous body at hand, we turn our attention to the physical principles governing the (thermo)mechanical behavior of continuum bodies. These principles are formulated in the form of *balance equations* of the (tensor) quantities under consideration. To give an overview of the most relevant equations, we formulate the general notation of a balance equation and subsequently use this framework to provide specific formulations for the relevant fields. For the case of singular surfaces, we refer the interested reader to the work by Šilhavý (1997). We assume in the following that there are no singular surfaces in the material body.

It is assumed that the change of a field quantity Λ is the sum of a production term <sup>Λ</sup> and a supply term <sup>Λ</sup> of Λ in the volume of the body ℬ and the flux <sup>Λ</sup> of Λ across the boundary *∂*ℬ, see Bertram (2012); Truesdell and Noll (2004) for more details. This relation can be expressed in integral form (Truesdell and Noll, 2004)

$$\frac{d}{dt} \int\_{\mathcal{B}\_t} \Lambda \, dv = \int\_{\mathcal{B}\_t} (p\_\Lambda + s\_\Lambda) \, dv + \int\_{\partial \mathcal{B}\_t} q\_\Lambda \cdot n \, da. \tag{2.19}$$

Here Λ, <sup>Λ</sup> and <sup>Λ</sup> denote tensor fields. <sup>Λ</sup> and <sup>Λ</sup> have the identical rank as Λ, whereas, for <sup>Λ</sup> the rank is increased by one. Using the divergence theorem in combination with Reynold's transport theorem enables to transform the integral form into a local form in regular points (Bertram, 2012)

$$\frac{\partial \Lambda}{\partial t} + \operatorname{div}(\Lambda v) = \operatorname{div}(q\_{\Lambda}) + p\_{\Lambda} + s\_{\Lambda}.\tag{2.20}$$

**Mass balance** The balance of mass emerges in local form by considering the density : ℬ × [0*,* ] → R. Thus, the balance of mass reads

$$
\dot{\rho} + \rho \operatorname{div}(v) = 0,\tag{2.21}
$$

as the production , supply and flux are zero for closed systems. For small strain deformations, it is typically assumed that the relation = <sup>0</sup> holds, i.e., the balance of mass is fulfilled for trivial reasons.

**Linear and angular momentum** For the linear momentum density , the supply term is given by the volume force density : ℬ ×[0*,* ] → R . The flux consists of the Cauchy stress tensor : ℬ×[0*,* ] → R × . Thus, the local form of the linear momentum balance reads (Šilhavý, 1997)

$$
\rho a = \operatorname{div}(\sigma) + b,\tag{2.22}
$$

where = ¨ denotes the acceleration, which, in the quasi-static setting considered in this work, is zero. The stress tensor describes the stress state in any material point. The symmetry of the stress tensor

$$
\sigma = \sigma^T \tag{2.23}
$$

is equivalent to the balance of angular momentum (Haupt, 2013).

**Energy and entropy** The sum of internal energy density : ℬ × [0*,* ] → R and kinetic energy density <sup>1</sup>*/*2 · is referred to as the total energy density. Thus, with the supply of internal heat sources : ℬ × [0*,* ] → R, the power of volume forces · , the negative heat flux − and the mechanical power · · , the second law of thermodynamics, i.e., the balance of total energy, reads in the local form

$$
\dot{e} + \frac{1}{2}\rho(v \cdot v) = b \cdot v + \omega + \operatorname{div}(\sigma^T \cdot v) - \operatorname{div}(q). \tag{2.24}
$$

The balance of internal energy results from subtracting the balance of linear momentum (2.22) multiplied by the velocity from equation (2.24) (Truesdell and Noll, 2004)

$$
\dot{e} = -\operatorname{div}(q) + w + \sigma : \dot{\varepsilon}.\tag{2.25}
$$

We use the entropy density : ℬ × [0*,* ] → R to formulate the local form of the entropy balance (Haupt, 2013)

$$\dot{s} = \text{div}(q\_s) + p\_s + s\_s,\tag{2.26}$$

where , and denote the flux, production and supply of entropy, respectively. According to the second law of thermodynamics, the entropy production is non-negative, i.e.,

$$p\_s \ge 0.\tag{2.27}$$

This conditions restricts the directions of thermodynamic processes (Coleman and Noll, 1974).

## **2.2 Deformation behavior of polycrystalline metals**

### **2.2.1 Crystal structure and deformation mechanisms**

With suitable tools to describe the motion of a continuum at hand, we focus on the material which is of primary interest in this work, i.e., the class of polycrystalline metals. They can be interpreted as an agglomerate of crystals or grains, see Fig.1.1. Each one of these crystals has a highly ordered and repetitive atomic lattice (Gottstein, 2013). Categorizing the different possible periodic arrangements of crystals results in the 14 different Bravais lattice types (Bravais, 1850). Due to its importance for industrial applications, we focus on the body-centered cubic (bcc) lattice, whose underlying unit cell is a cubic system with eight corner atoms and an additional center atom, i.e., in total two atoms per unit cell.

If a load is applied to the crystal lattice, at first it deforms elastically, i.e., the deformation is reversible, see Fig. 2.2b for a schematic drawing of an elastic deformation. When the load vanishes, the crystal lattice returns to its initial configuration, see Fig. 2.2a. It is experimentally observed that, if the applied load exceeds a material-specific limit, the structure deforms *plastically*, i.e., the deformation is irreversible. This process is depicted in Fig. 2.2, where the applied shear stress exceeds the limit <sup>∞</sup> resulting in a shift of atoms along a so-called glide plane by one inter-atomic distance.

**Figure 2.2:** Plastic deformation as the result of shifting the whole upper part at once by one atom. (a) = 0, (b)  *<* ∞, (c)  *>* ∞, (d) = 0.

Using the analogy of a crystal consisting of bubbles, Cottrell (1954) estimated the theoretical barrier for this deformation to be

$$
\tau\_{\infty} = \frac{\mu}{30},
\tag{2.28}
$$

where denotes the material's shear modulus. According to Hull and Bacon (2001) this theoretical limit for plastic deformation is, in general, several orders of magnitude *higher* than the experimentally observed yield strength.

An explanation for this striking difference was simultaneously proposed by Orowan (1934), Polanyi (1934) and Taylor (1934). Based on the presence of one-dimensional defects, so-called dislocations, their movement in the lattice, see Fig. 2.3a, is interpreted as the cause of plastic deformation. By breaking the atomic bonds of the inserted half plane with one of the neighboring planes and bonding again with the following edge, see Fig. 2.3b, the dislocation can move in the lattice. This movement facilitates the continuous plastic deformation of the crystal lattice and explains the reduced resistance against plastic deformation compared to equation (2.28). Peierls (1940) and Nabarro (1947) give a discussion of the resistance of a crystal lattice against dislocation movement.

**Figure 2.3:** Plastic deformation as a result of dislocation movement. (a) Initial configuration, (b) Moving dislocation, (c) Final configuration.

Dislocations generally do not move on planes with arbitrary normal, but on predefined systems depending on the crystal lattice (Hull and Bacon, 2001). The accumulation of dislocation movement on these systems is called slip. Those systems which permit dislocation movement are called slip systems (Hull and Bacon, 2001). For the body-centered cubic crystal there exist 12 primarily activated slip systems. Dislocation movement on other systems is only observed in some materials under certain temperature conditions (Courtney, 2005), which are not considered in this work. Within the crystal, a slip system is characterized by two orthogonal vectors, the slip plane direction and the normal .

For uniaxial tension applied to a crystal with slip systems , the driving force for dislocation movement, i.e., the shear stress, on the slip system follows Schmid's law (Gottstein, 2013)

$$
\tau^{\eta} = \sigma\_n \cos \phi^{\eta} \cos \lambda^{\eta}, \tag{2.29}
$$

where and denote the angles between the direction of applied stress and the slip plane normal and the slip direction, respectively.

For uniaxial tension, the shear stress becomes maximal for a slip system with = = 45<sup>∘</sup> .

### **2.2.2 Single crystal plasticity<sup>1</sup>**

To use the powerful computational homogenization methods outlined in Sec. 2.3, it is necessary to encapsulate the knowledge of the deformation behavior, see Sec. 2.2.1, into a suitable material model. In the case of polycrystalline materials, the crystal plasticity model offers a dedicated framework to describe the deformation behavior of polycrystalline materials accounting for the movements of dislocations. Based on the pioneering work by Asaro and Rice (1977), Asaro and Needleman (1985) and Peirce et al. (1983), many crystal plasticity models were developed, see Roters et al. (2010) for an overview.

Phenomenological crystal plasticity models use a small number of parameters to describe the observed deformation behavior of the crystallites. Nevertheless, they proved to be able to provide results which, for relevant features, agree well with experimental observations (Schäfer et al., 2019b; Natkowski et al., 2021b). For fatigue behavior, the deformations are typically small. Thus, we use a crystal plasticity model at small-strains, closely following Schäfer et al. (2019a). We assume the total strain tensor to be decomposed additively

$$
\varepsilon = \varepsilon\_e + \varepsilon\_p \tag{2.30}
$$

into an elastic and a plastic contribution, denoted by and , respectively. We refer to (Simo and Hughes, 2006, Cha. 1) for details. The Cauchy stress tensor , which is symmetric by conservation of angular momentum (Haupt, 2013), is assumed to be linear in the *elastic* strain,

<sup>1</sup> This subsection is based on excerpts from the publications by Kuhn et al. (2021; 2022)

i.e., Hooke's law

$$\sigma = \mathbb{C} : \varepsilon\_e \equiv \mathbb{C} : (\varepsilon - \varepsilon\_p), \tag{2.31}$$

involving the fourth-order stiffness tensor C, is assumed to hold. For elasto-viscoplasticity, the evolution of plastic strain is typically governed by a flow rule of the form

$$
\dot{\varepsilon}\_p = f(\sigma, z),
\tag{2.32}
$$

where denotes a vector of additional internal variables, see Chap. 2 in Simo and Hughes (2006).

In the framework of crystal plasticity, it is assumed that plastic deformation, i.e., ˙ ̸= 0, is caused by the movement of dislocations on the corresponding crystallographic slip systems. We consider volumepreserving slip mechanisms, i.e., conservative glide. An arbitrary slip system is activated when the resolved shear stress

$$
\tau^{\eta} = \sigma \cdot M^{\eta} \tag{2.33}
$$

exceeds a critical value . Here, denotes the symmetrized Schmid tensor of the slip system

$$M^{\eta} = \frac{1}{2} \left( m^{\eta} \otimes n^{\eta} + n^{\eta} \otimes m^{\eta} \right). \tag{2.34}$$

Following Bishop (1953) and Hutchinson (1976), the flow rule (2.32) is formulated as the superposition of crystallographic slip rates on the individual slip systems

$$\dot{\varepsilon}\_p = \sum\_{\eta=1}^{N\_S} \dot{\gamma}^{\eta} M^{\eta},\tag{2.35}$$

where ˙ denotes the plastic slip rate for the -th system. To determine the amount of plastic slip on a slip system , the theory of elastoviscoplasticity assumes the plastic slip rate ˙ to be a function of the resolved shear stress. The flow rule is designed in such a way that it captures, both, the Bauschinger effect and the ratcheting behavior, which are necessary for describing the cyclic behavior of metals accurately. The Bauschinger effect concerns a decreased resistance to plastic deformation also when a change of loading directions (e.g., from tension to compression) occurs. Ratcheting refers to a consistent accumulation of plastic slip under stress-driven cyclic loading conditions. Please note that the loading levels which lead to ratcheting under repeated loading will not induce a plastic flow of the material under static loading conditions, in general. For ratcheting and strain-driven loading, the mean stress will decrease for increasing number of cycles, see Farooq et al. (2020) and Cruzado et al. (2020) for recent discussions of ratcheting and the Bauschinger effect, respectively.

In this work, the flow rule proposed by Hutchinson (1976)

$$\dot{\gamma}^{\eta} = \dot{\gamma}\_0 \text{sgn}(\tau^{\eta} - \mathcal{X}\_b^{\eta}) \left| \frac{\tau^{\eta} - \mathcal{X}\_b^{\eta}}{\tau\_c^{\eta}} \right|^c \tag{2.36}$$

is used. The symbol denotes the critical resolved shear stress and refers to the stress exponent. To capture both the Bauschinger effect and ratcheting behavior, the original flow rule is augmented by the backstress , as proposed by Cailletaud (1992). Alternative approaches for modeling these effects are discussed by Harder (2001) with special focus on the coupling of the backstress evolution on different slip systems and the existence of the backstress tensor. Please note that the vector of internal variables in equation (2.32) collects the different backstresses .

To close the model, the flow rule needs to be complemented by an evolution equation for the backstress . The Armstrong-Frederick (AF) model (Frederick and Armstrong, 2007) is a nonlinear extension of the model by Prager (1949) and involves a recall term. The model describes

the evolution of the backstress by the equation

$$\dot{\mathcal{X}}\_b^\eta = A \,\dot{\gamma}^\eta - B \,\mathcal{X}\_b^\eta \, |\dot{\gamma}^\eta|\,, \tag{2.37}$$

where and are material parameters with dimensions MPa and 1, respectively. These parameters represent the direct hardening as formulated in the Prager model and the recovery modulus present in the extended model, respectively. Bari and Hassan (2000) showed that the AF model may overestimate ratcheting encountered in experiments. Furthermore, they attribute this shortcoming to the constant ratcheting associated to the model.

Extending a formulation by Chaboche, which is based on a decomposition of the AF model, Ohno and Wang (1993) proposed a kinematic hardening model that reads, in the context of crystal plasticity,

$$\dot{\mathcal{X}}\_b^\eta = A \dot{\gamma}^\eta - B \left(\frac{|\mathcal{X}\_b^\eta|}{A/B}\right)^M \mathcal{X}\_b^\eta \left|\dot{\gamma}^\eta\right|. \tag{2.38}$$

In the Ohno-Wang (OW) model, the dynamic recovery term is modified by a power law with exponent , which controls the degree of non-linearity, to improve the prediction of ratcheting and mean stress relaxation.

According to Schäfer et al. (2019a;b) and Natkowski et al. (2021a;b) the above model captures relevant aspects of fatigue. Thus, in the following it serves as our model for the deformation behavior of polycrystalline microstructure.

## **2.3 Homogenization methods**

### **2.3.1 Mechanical properties of microstructured materials**

The deformation mechanisms on the microscopic scale, e.g., outlined in Sec. 2.2.1 for polycrystalline metals, influence the overall macroscopic behavior of the material. Thus, within many applied materials, we face different scales when investigating mechanical properties. Let us consider the component shown in Fig. 2.4 as an example. On the component, or macroscopic, scale, the metallic material is assumed to be homogeneous. On the microscopic scale, see Fig. 2.4 for a example of the components microstructure, we observe that this not the case. Owing to the spatial distribution of crystallographic orientations and the direction-dependent mechanical properties of crystallites, the microstructure is *heterogeneous*. These microstructural heterogeneities influence the macroscopic mechanical properties and can be crucial for failure of the component, see François (1989) and Krupp (2007) for an overview in the context of fatigue failure.

Deriving macroscopic properties based on information on the microscale is called *homogenization* whereas localization refers to deducing microstructural behavior from macroscopic quantities, e.g., strains and stresses. To determine the effective properties of a microstructure, we have to consider a domain which is *characteristic* for the whole microstructure. This representative volume element (RVE) (Hill, 1952; 1963) needs to be sufficiently large to capture all relevant statistics of the microstructure, while being sufficiently small compared to the overall size of the specimen. We show this concept of scale separation schematically in Fig. 2.4. The component scale *ℓ*macro is much larger than the microstructural scale *ℓ*micro, which is larger than scale of one constituent, *ℓ*constituent, e.g., the size of one grain in polycrystalline metals.

**Figure 2.4:** Schematic illustration of different lengths on the micro and macroscopic scale.

Concluding that the mechanical behavior of the RVE is representative for the whole material allows determining effective macroscopic fields, e.g., strains ¯ and stresses ¯, from averaging the fluctuating micro stress and strain fields, () and () over the RVE volume , i.e.,

$$
\bar{\sigma} = \frac{1}{V} \int\_{V} \sigma(x) \, dV \tag{2.39}
$$

$$
\bar{\varepsilon} = \frac{1}{V} \int\_{V} \varepsilon(x) \,dV. \tag{2.40}
$$

To keep this section focused on the main aspects in this work, we discuss analytical and computational homogenization methods with their respective requirements as well as a advantages and disadvantages in the following. For more details regarding micromechanics and homogenization methods we refer to Nemat-Nasser and Hori (1993), Torquato (2002) and Gross and Seelig (2017).

#### **2.3.2 Analytical approaches**

The works of Voigt (1889) and Reuss (1929) describe the first analytical approaches to homogenization. Both authors approximate the effective elastic moduli of a polycrystal by volume averaging the elastic constants of single phases.

Eshelby (1957) provided an analytical solution for the strain and stress fields of an eigenstrained ellipsoidal inclusion in an infinite homogeneous matrix, for which the strains and stresses inside the inclusion are homogeneous. Based on this solution, the mean-field theory provides many different homogenization approaches. These methods have successfully been applied to determine effective properties of fiber reinforced composites (Schemmann et al., 2018; Kehrer et al., 2020) and polycrystalline materials (Rieger and Böhlke, 2015). For the linear elastic effective behavior of polycrystals, the self-consistent method, proposed by Hershey (1954) and Kröner (1958), is used frequently. There exist many different extensions to apply this method for non-linear behavior as well, e.g., (Hill, 1965; Molinari et al., 1987). In their article, Lebensohn et al. (2007) give an overview over various material classes they study using the self-consistent approach. Applying a variational framework, Hashin and Shtrikman developed bounds on the effective properties of heterogeneous microstructures (Hashin and Shtrikman, 1962a;b). These bounds are used in a variety of applications, e.g., Fernández and Böhlke (2019) applied the Hashin-Shtrikman bounds to the effective elastic properties of textured polycrystals. We refer the interested reader to Chap. 8 in Gross and Seelig (2017) for a more detailed overview over analytical homogenization methods.

Although some restrictions apply to these techniques, e.g., in terms of the possible inclusion geometries considered in Eshelbys solution, they are capable of predicting the effective mechanical properties of various material classes. Additionally, they are efficient both in terms of

computing time and memory requirements, making them attractive for industrial applications.

#### **2.3.3 Computational homogenization**

To circumvent the simplifying assumptions of analytical approaches, e.g., regarding the inclusion geometry, computational homogenization approaches use spatially resolved microstructures as computational cells. In addition to the effective behavior, numerically solving the so-called cell problem on the resolved microstructure, delivers the full fields of stresses and strains, which emerge. Furthermore, it is possible to compute the resulting plastic deformations developing within the microstructure on an intracrystalline resolution which can be used as input, e.g., to predict crack initiation and propagation in polycrystalline metals (Schäfer et al., 2019a;b; Natkowski et al., 2021a;b).

To derive the emerging fields and thus the effective properties of a microstructure, the finite element method (FEM) might be used to solve the discretized weak form of the governing equations (Matouš et al., 2017). This method predicts the effective properties, e.g., of composites (Nakamura and Suresh, 1993; Kouznetsova et al., 2001; Borbely et al., 2001), microstructures with pores (Zhuang et al., 2015) and polycrystals (Miehe et al., 1999; Shenoy et al., 2008; Vajragupta et al., 2014).

In this work we use a method based on the work by Moulinec and Suquet (1994; 1998) to numerically determine the effective properties of polycrystalline microstructures. This approach is based on the reformulation of the cell problem into a Lippmann-Schwinger equation and uses the fast Fourier transforms (FFT) for its computational resolution. Due to its low memory footprint, the efficient FFT implementations and the straightforward realization, the method became increasingly popular over the past decade. This fueled the development of new algorithms, see Schneider (2021) for a recent review on supported discretization methods and available solution schemes. Applications of this method

include, but are not limited to, damage and fracture mechanics (Chen et al., 2019; Ernesti et al., 2020), homogenization of piezoelectricity (Brenner, 2009; Göküzüm et al., 2019), computing the deformation behavior of polycrystalline materials (Lebensohn and Rollett, 2020; Rovinelli et al., 2020) and composites (Müller et al., 2015; Görthofer et al., 2020). We rely on FeelMath (Fraunhofer ITWM, 2021), a commercial FFT based micromechanics solver, to determine a solution of our cell problem. Fig. 2.5 shows one example for the resulting von Mises stress field in a polycrystalline microstructure computed by FeelMath. Nevertheless we still need to supply computational cells, i.e., RVEs of the considered polycrystalline microstructure. Additionally we need to identify parameters for the material model described in Sec. 2.2.2 to describe the deformation behavior on a microscopic scale accurately.

**Figure 2.5:** Example for computational homogenization. (a) Spatially resolved polycrystalline microstructure, (b) Resulting von Mises stress field at 0*.*1% macroscopic strain.

## **2.4 Prerequisites for computational homogenization of polycrystalline metals**

### **2.4.1 Polycrystalline computational cells for micromechanical simulations**

The morphology of polycrystalline microstructures can be determined by dedicated experimental methods, like serial sectioning (Spowart et al., 2003; Kubis et al., 2004) and FIB-SEM (Bansal et al., 2006; Groeber et al., 2006; Zaefferer et al., 2008)). By combining these with Electron Backscatter Diffraction (EBSD) methods (Adams and Olson, 1998; Korte et al., 2011) the local orientation of each grain can be resolved. Fig. 1.1 and Fig. 2.4 show two examples of images obtained by EBSD.

Directly using these images as an input for computational homogenization might seem like the most natural approach. However, the experimental procedures are quite time and cost intensive (Bargmann et al., 2018) and, in addtion, determining a *representative* image for the whole microstructure can be challenging. Due to size limitations of one image, it might not be large enough to capture statistical fluctuations of the relevant microstructure features. It might even distort these. Using the (generally) non-periodic experimental images can lead to large required computational cells, e.g., in terms of grain number within the image, to accurately predict the effective behavior (Kanit et al., 2003; Yang et al., 2019; Schneider, 2021).

As an alternative, synthetically generated computational cells gained popularity, see Fig. 2.6 for a polycrystalline example. These synthetic microstructures aim at reproducing statistical features of the real microstructure, e.g., grain size and aspect ratio distribution (Gillner and Münstermann, 2017; Prasad et al., 2019; Henrich et al., 2020). They offer the benefit of direct control over the microstructural parameters, yielding the possibility to tailor the microstructure to certain loading

conditions and to capture the statistical influence of fluctuations in the microstructural descriptors (Tallman et al., 2020; Tran and Wildey, 2021).

**Figure 2.6:** Example of a polycrystalline computational cell.

There exists a variety of methods to generate synthetic microstructures, each tailored to the material class under investigation. In the following we focus on polycrystalline metals and refer the interested reader to the review article by Bargmann et al. (2018) for a general overview.

One possible approach is simulating the solidification and grain growth occurring during the manufacturing process of a metal. Examples include techniques based on Vertex methods (Kawasaki et al., 1989; Weygand et al., 1999), cellular automata (Janssens, 2010; Saluja et al., 2012; Kühlbach et al., 2016) and Monte Carlo methods (Anderson et al., 1989; Radakrishnan and Zacharia, 1995). Using the phase field method (Krill III and Chen, 2002; Nestler, 2005; Hoetzer et al., 2016) allows considering physical principles when simulating the forming process of metals. This comes at the cost of increased computing time and having to identify physically reasonable model parameters. Examples for polycrystalline

materials using the phase field approach can be found in the work by Graf et al. (2021a;b), simulating the formation of martensite from a prior austenite phase for an industrial application.

Describing the grains as geometric objects dispenses with the computational demanding simulation of solidification and grain growth. These objects form a tessellation of the computational domain, i.e., there are no voids and the objects do not overlap. Assuming the grains to be convex polyhedrons, Voronoi tessellations (Aurenhammer, 1991) can be used to generate representations of polycrystalline microstructures. This approach uses a distinct set of points, called seeds, in the computational domain to define a cell by all the points lying closer to one seed than to any other seed, see Fig. 2.7 for a two dimensional example.

**Figure 2.7:** Schematic drawing of a Voronoi tessellation process, following Aurenhammer (1991). (a) Seeds in the computational domain, (b) Resulting tessellation.

For seed positions following a Poisson process, a so-called Poisson-Voronoi tessellation, the generated microstructures often do not capture relevant properties of the real microstructure, e.g., grain size distribution or number of neighboring grains (Döbrich et al., 2004; Luther and Könke, 2009). Choosing the seeds based on a precomputed sphere packing increases the synthetic microstructures realism (Fritzen et al., 2009; Quey et al., 2011). Laguerre tessellations extend the Voronoi tessellations by one additional degree of freedom per cell. This so-called tessellation weight controls the size of each cell, enabling the generation of polycrystalline microstructures with prescribed volume fractions (Bourne et al., 2020). For generating synthetic microstructures with non-convex geometric models, several approaches are available, e.g., fitting experimental data (Alpers et al., 2015; Teferra and Rowenhorst, 2018) or optimization procedures to match the shape distributions (Groeber and Jackson, 2014; Prasad et al., 2019; Henrich et al., 2020).

The heterogeneity of a polycrystalline microstructure is not only influenced by the shape of its grains, but especially by the distribution of crystallographic orientations (Adams and Olson, 1998; Kocks et al., 2000). Thus, the geometric modeling approaches require an additional step, the assignment of crystallographic orientations to each grain model. For polycrystalline materials it is important to distinguish between two cases: The uniformly distributed orientations and the textured case where certain orientations are more likely to be encountered than others. The former case results in mechanically isotropic effective properties (Krawietz, 1999; Bertram et al., 2000; Böhlke and Bertram, 2001), whereas the elastic and elasto-plastic behavior in the latter case might, in general, be anisotropic. Besides randomly assigning orientations, e.g., see (Gillner and Münstermann, 2017; Tu et al., 2019), some dedicated algorithms exist. Examples include approaches discretizing the one-point correlation function, i.e., the crystallite orientation distribution function (CODF), which describes the orientation state of a polycrystal (Tóth and Van Houtte, 1970; Melchior and Delannay, 2006). Following their respective discrete CODF values, crystallographic orientations are selected and assigned to cells in the synthetic microstructure (Chunlei et al., 2000; Deka et al., 2006). Building upon this discretization, Eisenlohr and Roters (2008) use a combination of stochastic and deterministic approaches

to increase sampling efficiency. Biswas et al. (2020) and Prasad et al. (2019) iteratively adjust the discretization of the CODF to improve the sampling results, whereas Liu et al. (2020a) adjust the parameters of the microstructure generator. Quey and Renversade (2018) match the special case of a uniform orientation distribution using uniformly distributed unit quaternions. Approximating the CODF by superposing so-called texture components (Helming, 1996; Wassermann and Grewen, 2013) allows sampling these components and assigning orientations accordingly (Hielscher and Schaeben, 2008).

### **2.4.2 Determination of single crystal plasticity parameters**

To describe the deformation behavior of polycrystalline metals on a microstructural level, we use the crystal plasticity model outlined in Sec. 2.2. However, this model still needs to be equipped with suitable parameters. If the model parameters cannot be derived using physical principles, we will need to determine the parameters so that the model output matches experimental observations. Thus, we need to minimize a function which describes the difference between experiment and simulation.

For crystalline metals, ultrasonic testing on a single crystal (Huntington, 1947), microtensile testing (Gianola and Eberl, 2009) as well as experiments using milled micropillars (Cruzado et al., 2015) and nanoindentation tests (Chakraborty and Eisenlohr, 2017; Engels et al., 2019; Liu et al., 2020b) give insight into the deformation behavior on a microscopic scale. Comparing the experimental observations with simulations on single crystals, in combination with suitable optimization methods, allows determining the crystal plasticity parameters (Chakraborty and Eisenlohr, 2017; Engels et al., 2019).

Using data from *polycrystalline* experiments allows reducing the experimental effort associated with single crystal tests. Comparing the macroscopic data with the effective behavior, determined from simulations

on polycrystalline RVEs, allows calibrating crystal plasticity parameters, see Herrera-Solaz et al. (2014). Schäfer et al. (2019a) define a cost function based on the difference between experimental and simulated stress-strain hysteresis, leading to an optimization problem where the objective function is expensive to compute and gradient information is not available. Because automatic differentiation (Griewank and Walther, 2008) might be difficult to integrate into existing homogenization tools, derivative free approaches are frequently used in this setting, see Rios and Sahinidis (2013) for an overview. However, care has to be taken when considering industrial applications, as determining suitable parameters are a regular task (Schäfer et al., 2019a;b; Natkowski et al., 2021b; Arnaudov et al., 2020). Indeed, as evaluating the cost function requires a full-field simulation it is desirable to minimize the number of function evaluations.

**Chapter 3**

# **Fast methods for computing centroidal Laguerre tessellations for prescribed volume fractions with applications to microstructure generation of polycrystalline materials <sup>1</sup>**

## **3.1 Introduction**

In a recent article, Bourne et al. (2020) used a convex variational principle for generating Laguerre tessellations modeling polycrystalline microstructures with prescribed volume fractions per grain, see Fig. 3.1 for an overview. Their insights led to the development of *fast* techniques for generating polycrystalline microstructures which are orders of magnitude faster than other approaches in the literature, for instance based on derivative-free optimization (Quey and Renversade, 2018), statistical methods (Lautensack, 2008) or crossentropy approaches (Petrich et al., 2019), and more accurate than

<sup>1</sup> This chapter is based on the publication by Kuhn et al. (2020) whereas minor typographical and formatting changes have been made for cohesion of the manuscript and the readers convenience

heuristic approaches (Lyckegaard et al., 2011).

In this chapter we improve upon the proposed method by Bourne et al. (2020), as they essentially rely upon black-box optimization methods. They use the limited-memory Broyden-Fletcher-Goldfarb-Shanno (BFGS) method (Nocedal, 1980) with backtracking for solving the convex optimization problem. The backtracking is used to avoid generating empty cells during the iterative process, and is also used for the corresponding Newton-type methods (Kitagawa et al., 2019).

In section 3.2, we set up our notation for computing periodic Laguerre tessellations and the resulting convex optimization problem. We discuss pertinent solvers in section 3.3. Of particular interest here are modern non-monotone gradient-type solvers, i.e., the methods of Malitsky and Mishchenko (2019) and Barzilai and Borwein (1988), whose application to optimal-transport problems appears to be novel. Compared to more traditional solvers of (Quasi-)Newton type, they combine low cost per iteration with non-monotone convergence behavior, which turns out to be beneficial for the *overall* run-time. Indeed, non-monotone methods have gained popularity in recent years (Nesterov, 2004), as the user is mostly interested in minimizing the time to solution, and ensuring a monotone convergence of the algorithm in question is typically of secondary importance.

For computing centroidal Laguerre tessellations, we review the classical Lloyd's algorithm in subsection 3.2.2. Subsequently, we apply Anderson acceleration to Lloyd's algorithm. Anderson acceleration (Anderson, 1965) is a dedicated technique for improving the convergence behavior of general fixed-point iterations. Its mathematical analysis is more recent, establishing it as a particular Quasi-Newton method of multi-secant type (Fang and Saad, 2009) and establishing the improved convergence behavior also rigorously (Toth and Kelley, 2015; Evans et al., 2020).

In section 3.4, we study the numerical effectiveness of the proposed algorithms for selected problems of increasing complexity. We study the dependence of the algorithms on the number of grains and the prescribed grain-size distribution, with and without centering. Furthermore, we exhibit their capabilities for polycrystalline materials of industrial relevance.

**Figure 3.1:** General flowchart for computing Laguerre-type models of polycrystalline microstructures, Bourne et al. (2020).

## **3.2 Periodic Laguerre tessellations and centering**

#### **3.2.1 Periodic Laguerre tessellations**

Let be a rectangular domain in R , which we fix as = [0*,* 1) × [0*,* 2) × *. . .* × [0*,* ) for positive lengths . shall serve as our unit cell of interest, and we furnish it with periodic boundary conditions. We

denote the corresponding periodic distance function by : × → R≥0. For a sequence (1*,* 2*, . . . ,* ) ∈ of distinct points lying inside and a given set of real weights ∈ R , we denote by 1*,* 2*, . . .*  the periodic Laguerre tessellation (Aurenhammer, 1991) of with seeds (1*,* 2*, . . . ,* ) and weights (1*,* 2*, . . . ,* ). Each cell is a subset of defined via

$$D\_i = \left\{ x \in Y \, \middle| \, d(x, x\_i)^2 - w\_i \le d(x, x\_j)^2 - w\_j \quad \text{for} \quad j \in \{1, 2, \dots, K\} \right\}. \tag{3.1}$$

The interiors of the Laguerre cells are mutually disjoint, and the Laguerre cells cover the original volume . Clearly, adding a constant to each single weight does not change the associated Laguerre cells. Notice that the seed does not need to lie in the corresponding Laguerre cell . Also, the Laguerre cell may be empty. This becomes clear when imagining the weight as a "price". Loosely speaking, if the price is excessive, the turnover, which is related to the volume of the cell, will be zero.

For the special case of identical weights, <sup>1</sup> = <sup>2</sup> = *. . .* = , a Laguerre tessellation is called a Voronoi tessellation, see Aurenhammer (1991). For this special case, the seeds are always contained in the corresponding cell. However, the kinds of tessellations which may be generated as Voronoi tessellations are rather limited, for instance in terms of the achievable volume-fraction distributions, see section 2.1.1 in Quey and Renversade (2018). In contrast, Laguerre tessellations may generate any normal tessellation with convex cells in R for ≥ 3, see Lautensack and Zuyev (2008). In this context, normal means that if the seeds are in general position, each -face of the tessellation arises as the intersection of exactly − + 1 cells for = 0*,* 1*, . . . ,*  − 1.

For this article, Laguerre tessellations serve as simple geometrical models for polycrystalline microstructures, e.g., of metallic materials, alloys, ceramics or polycrystalline ice. We refer to section 2.1 in Bargmann et al. (2018) for a recent overview.

For a given sequence (1*,* 2*, . . . ,* ) ∈ of distinct points of lying inside and a given set of volume fractions (1*,* 2*, . . . ,* ) ∈ R , each of which is non-negative, and which sum to unity, there is a corresponding set of weights (1*,* 2*, . . . ,* ) ∈ R , s.t. the corresponding Laguerre cells realize the prescribed volume fractions, i.e.,

$$\frac{|D\_i|}{L\_1 L\_2 \cdots L\_d} = \phi\_i,\tag{3.2}$$

see Aurenhammer et al. (1998), where | | denotes the volume of cell . The proof of this result is interesting in its own right, as it immediately gives rise to a computational technique for computing the associated weight vector . As noted earlier, the component-wise addition of a constant to the weight vector does not change the associated Laguerre cells. Thus, we define

$$W\_K = \left\{ w \in \mathbb{R}^K \left| \sum\_{i=1}^K w\_i = 0 \right. \right\}, \tag{3.3}$$

the ( − 1)-dimensional subspace of R with vanishing mean. For a fixed sequence (1*,* 2*, . . . ,* ) ∈ of distinct points lying inside and a given set of volume fractions (1*,* 2*, . . . ,* ) ∈ R , define the function

$$\begin{aligned} g: W\_K &\to \mathbb{R}, \\ w &\mapsto \sum\_{i=1}^K \left( \phi\_i \cdot L\_1 L\_2 \cdots L\_d - |D\_i| \right) w\_i + \sum\_{i=1}^K \int\_{D\_i} d(x, x\_i)^2 \, dx, \end{aligned} \tag{3.4}$$

where denotes the -th associated Laguerre cell (3.1). Aurenhammer-Hoffmann-Aronov (Aurenhammer et al., 1998) have shown that the function (3.4) is a concave function which admits maximizers. Furthermore, each critical point (which is a maximum by concavity of ) satisfies the first order optimality condition

$$0 = Dg(w)[v]$$

in terms of the directional derivative of at ∈ in direction ∈ . More explicitly,

$$\left.Dg(w)[v]\right| \equiv \frac{d}{ds}g(w+sv)\Big|\_{s=0} = \sum\_{i=1}^{K} \left(\phi\_i \cdot L\_1 L\_2 \cdots L\_d - |D\_i|\right) v\_i,\tag{3.5}$$

i.e., any critical point of satisfies equation (3.2). Furthermore, it can be shown that is even strictly concave on the set

$$\left\{ w \in W\_K \, | \, |D\_i| > 0 \quad \text{for all} \quad i = 1, 2, \dots, K \right\},\tag{3.6}$$

see Kitagawa et al. (2019). It is clear that is not strictly concave outside of the set (3.6). Indeed, suppose that we have ∈ , s.t. the -the Laguerre cell is empty. Then, the Laguerre cells remain unchanged if we increase the -th weight . In particular, the function is constant along the ray

$$w + s \left[ e\_i - \frac{1}{K} (1, 1, \dots, 1) \right], \quad s \ge 0, 1$$

where is the -th unit vector in R . For completeness, notice that the -th component of the gradient of on computes as

$$\nabla\_i g(w) = \phi\_i \cdot L\_1 L\_2 \cdots L\_d - |D\_i|,\tag{3.7}$$

see (3.5). Notice that ∇() ∈ automatically, i.e., its mean vanishes. The concavity of (or, equivalently, the convexity of −) permits using established convex optimization techniques as in Nocedal and Wright (1999) to be applied, see section 3.3 for details. To finish this paragraph, we report on a result in Bourne and Roper (2015). The function (3.4) is

even twice differentiable on with second derivative

$$D^2g(w)[v,v] = v^T H(w)v, \quad w, v \in W\_K,$$

in terms of the Hessian matrix () ∈ R × with entries

$$H\_{ij}(w) = \frac{A\_{ij}}{2\,d(x\_i, x\_j)}\tag{3.8}$$

for ̸= , where denotes the area of the face between the -th and -th cell (which is formally set to zero if the cells do not share a face), and = − ∑︀ ̸= for = 1*,* 2*, . . . ,* .

As is concave and twice differentiable, the Hessian matrix () is negative semidefinite. However, () does not have full rank, as all rows and columns have zero mean, leading to a zero eigenvalue. In particular, it is not invertible. As a countermeasure, we introduce the matrix

$$
\tilde{H}(w) = H + \frac{\text{tr}(H(w))}{K^2} \underline{1}\,\underline{1}^T \tag{3.9}
$$

with the vector 1 = (1*,* 1*, . . . ,* 1) ∈ R of all ones. The latter modification makes sure that

1. If the system

$$
\tilde{H}(w)w = v
$$

is solved for Δ ∈ R with ∈ , then Δ ∈ automatically.

2. If the eigenvalues of (), considered as an operator on , are contained in [−*,* +], the eigenvalues of ˜ () are contained in [−*,* +] as well. Furthermore, if | | *>* 0 for all = 1*,* 2*, . . . ,* , the function is even strictly convex, and, thus, the eigenvalues of () on are strictly negative, i.e., <sup>+</sup> *<* 0. As a consequence, also the eigenvalues of ˜ () are strictly negative, ensuring invertibility of ˜ ().

These two properties are useful when implementing Newton's method, see section 3.3.4. Item 1 means that we may essentially forget about the space , and work in R without second thoughts. The second item makes sure that our modification does not lead to an ill-conditioning of the linear system determining (3.24) the Newton increment. Indeed, for any  *>* 0, the matrix

$$H - \lambda \underline{1} \underline{1}^T$$

is negative semi-definite. If is not chosen carefully, it may deteriorate the conditioning of the matrix in question. We know that

$$\frac{\text{tr}(H(w))}{K} \in [\lambda\_-, \lambda\_+],$$

as it corresponds to the mean of the eigenvalues. Furthermore, the unit vector 1*/* √ 1 spans the kernel of () (at least on the set where is strictly concave), and 1*/* 1 1 is the orthogonal projector onto this kernel. Combining these observations leads to item 2.

#### **3.2.2 Generating centroidal Laguerre tessellations**

For a tessellation {} of denote by ∈ a centroid of , i.e., a minimizer of

$$\int\_{D\_i} d(c, x)^2 \, dx \longrightarrow \min\_{c \in Y} \dots$$

The periodic boundary conditions may give rise to special situations if the cells "wrap around ". As a simple example, for the single cell , *any* ∈ will be a centroid.

However, if the Laguerre cell is sufficiently small, i.e., it may be displaced in to a configuration not touching the boundary *∂* ⊆ R , the classical representation as the center of mass

$$c\_i = \frac{1}{|D\_i|} \int\_{D\_i} x \, dx \tag{3.10}$$

(relative to a suitable local coordinate system) may be used. We shall rely upon the latter expression for practical purposes.

A Laguerre tessellation is called centroidal (Brieden and Gritzmann, 2012) if each seed is a centroid of its corresponding cell. Although, when looking at solidification of polycrystalline materials, there does not seem to be any requirement that the nucleation site corresponds to the grain centroid, it has been observed, see Bourne et al. (2020), that centroidal Laguerre tessellations exhibit a better shape regularity than general Laguerre tessellations.

Suppose a set of (mutually different) seeds (1*,* 2*, . . . ,* ) ∈ and volume fractions ∈ R with *>* 0 and ∑︀ =1 = 1 are given. By the reasoning of the previous section, we find a set of weights ∈ R whose associated Laguerre tessellation with seeds (1*,* 2*, . . . ,* ) has cells with volume fractions . Denote by (1*,* 2*, . . . ,* ) ∈ the centroids of the respective Laguerre cells. A simple method for iteratively approaching a centroidal Laguerre tessellation with volume fractions as above is to use the centroids as new seeds. A pseudo-code of this classical Lloyd-type algorithm (Bourne and Roper, 2015) is shown in Alg. 1.

#### **Algorithm 1** Lloyd-type centering (*,* maxit*,* tol)

1: Determine initial guess ( 0 1 *,* <sup>0</sup> 2 *, . . . ,* <sup>0</sup> ) of seeds


```
6: compute centroids (

                            1
                             , 
                               2
                                , . . . , 
                                       ) by (3.10)
```

```
7: (
        +1
        1
            , +1
              2
                  , . . . , +1
                          ) ← (

                                    1
                                     , 
                                       2
                                        , . . . , 
                                               )
```

```
8:  ←  + 1
```

Several remarks are in order:


$$\text{residual} = \sqrt{\frac{\sum\_{i=1}^{K} \phi\_i \, d(x\_i^k, c\_i^k)^2}{(\max L\_i)^2}},\tag{3.11}$$

i.e., the 2 -norm of the difference between seeds and centroids, weighted by volume-fraction and normalized by the maximum edge length of the unit cell. A typical value for tol is 10<sup>−</sup><sup>1</sup> .

4. The Lloyd-type algorithm presented can be shown to converge by using a monotonicity argument (Bourne and Roper, 2015).

The problem of obtaining a centroidal Laguerre tessellation may also be written as an energy-minimization problem, see Liu et al. (2016). However, the resulting optimization problem admits a multitude of minimizers, i.e., we face a non-convex optimization problem, in general, and it is difficult to design powerful computational strategies. For instance, there are examples where the BFGS method does not converge, even with exact line search, see Dai (2012).

As an alternative, we investigate Anderson acceleration (Anderson, 1965), which is a technique for accelerating general fixed point iterations

$$x^{k+1} = F(x^k),$$

where, in our context, ∈ and the mapping

$$F: Y^K \to Y^K \tag{3.12}$$

is given by computing the centroids of the Laguerre tessellation, see Alg. 1.

More precisely, the Anderson-accelerated version of the fixed point iteration (3.12) with depth  *>* 0 stores the last iterates

$$x^k, x^{k-1}, \dots, x^{k-m+1}$$

and their images w.r.t.

$$F(x^k), F(x^{k-1}), \dots, F(x^{k-m+1}),$$

computes ∈ R as the minimizer of the problem

$$\left\| \sum\_{i=1}^{m} \alpha\_i (F(x^{k+1-i}) - x^{k+1-i}) \right\|^2 \longrightarrow \min\_{\sum\_{i=1}^{m} \alpha\_i = 1} \tag{3.13}$$

and proposes the next iterate as

$$x^{k+1} = \sum\_{i=1}^{m} \alpha\_i F(x^{k+1-i}).\tag{3.14}$$

Several remarks are in order.


there are countably many vectors ∈ R with

$$y\_i = z\_i + \xi\_i \mod \quad L\_i \quad (i = 1, 2, \dots, d). \tag{3.15}$$

Indeed, for any such , and any -tuple (1*,* 2*, . . . ,* ) of integers, we may form the vector

$$
\tilde{\xi} = \xi + \sum\_{i=1}^{d} k\_i L\_i e\_i,\tag{3.16}
$$

which also satisfies (3.15) with replaced by ˜ . In turn, any two vectors and ˜ satisfying equation (3.15) are related by (3.16) by a suitable tuple (1*,* 2*, . . . ,* ) of integers. Thus, we define − as the vector satisfying (3.15) with minimal length (which is welldefined if and are sufficiently close in the periodic distance). Using this definition -times clarifies how the difference vectors ( +1− ) − +1− in (3.13) have to be understood.


For degenerate problems, safeguarding strategies are necessary for improving convergence (Fu et al., 2020).

7. The Anderson-accelerated Lloyd-type centering generalizes the relaxed (or damped) Lloyd method for = 1 with variable relaxation.

The resulting algorithm is summarized in Alg. 2. For = 0, we (formally) recover the original algorithm Alg. 1. For our numerical experiments, we set = 4.

#### **Algorithm 2** Anderson-accelerated Lloyd-type centering (*, ,* maxit*,* tol)

```
1: Determine initial guess (
                            0
                            1
                             , 0
                               2
                                , . . . , 0
                                       ) of seeds
2:  ← 0
3: residual ← +∞
4: while  < maxit and residual > tol do
5: compute Laguerre tessellation with seeds (

                                                 1
                                                  , 
                                                     2
                                                      , . . . , 
                                                            )
      for volume fractions 
6: compute centroids (

                           1
                            , 
                              2
                               , . . . , 
                                     ) by (3.10)
7: update the linear system for solving (3.13)
8: determine the vector 
                             by solving (3.13)
9: update 
               +1
                ←
                      ∑︀
                        =1 (
                                 +1−
                                       )
10:  ←  + 1
11: update residual, see equation (3.11)
12: end while
13: return ,residual
```
## **3.3 Optimization methods**

In this section, we describe pertinent optimization algorithms used for solving the concave optimization problem (3.4)

$$g(w) \longrightarrow \max\_{w \in W\_K}$$

for the unknown weight vector, and given seeds and volume fractions (which we assume to be fixed throughout this entire section). To bring it into a more convenient form, we shall equivalently rewrite the latter problem as a minimization problem via

$$f(w) \longrightarrow \min\_{w \in W\_K} \tag{3.17}$$

with () = −(), i.e., in terms of the convex function . In particular, we shall use descent methods instead of ascent algorithms, for instance.

#### **3.3.1 The gradient descent method**

Recall that the gradient of the function (3.17) is given, see equation (3.7), component-wise by

$$\nabla\_i f(w) = |D\_i| - \phi\_i \cdot L\_1 L\_2 \cdots L\_d \tag{3.18}$$

Thus, a simple gradient descent method

$$w^{n+1} = w^n - \alpha\_n \nabla f(w^n),\tag{3.19}$$

involving a positive sequence of step sizes , may be used for solving (3.17), as pioneered by Aurenhammer et al. (1998). If the (local) Lipschitz constant of the gradient of is known to be , the iterative scheme (3.19) can be shown to converge for the choice = 1*/*, see Nesterov (2004). However, for convex , the convergence rate is only logarithmic, in general. If, in addition, is -strongly convex in the vicinity of the minimizer, the convergence rate improves to be linear, and the choice = 2*/*( + ) optimizes the convergence rate (estimate). Notice that the required iteration count is proportional to */*.

Unfortunately, for the problem (3.17) both the Lipschitz constant and the strong convexity constant are not easily available. Of course, both may be obtained as the minimum and maximum eigenvalue of the Hessian (3.8). However, even computing the entries of the Hessian matrix (not even accounting for extracting the extremal eigenvalues) leads to significant computational overhead.

We are only interested in the qualitative performance of the gradient scheme (3.19) in order to comparing it asymptotically to more elaborate minimization schemes. Thus, taking into account the necessary dimension of the step-size (the weights have dimensions length<sup>2</sup> , the gradient has dimension ), we simply use

$$
\alpha\_n = 0.1 \left( L\_1 L\_2 \cdots L\_d \right)^{\frac{2}{d} - 1}. \tag{3.20}
$$

Aurenhammer et al. (1998) propose using a step-size of Polyak-type (Polyak, 1987). However, the latter requires a good estimate of the optimal objective-value (which is not available, in general) to be competitive. Thus, we stick to the simple choice (3.20).

#### **3.3.2 The Malitsky-Mishchenko solver**

Malitsky and Mishchenko (2019) designed an adaptive strategy for choosing the step size in gradient methods which is globally convergent for general convex continuously differentiable functions (not necessarily with globally Lipschitz-continuous gradient) that automatically exploits strong convexity. The method is non-monotone, but significantly outperforms gradient descent methods with fixed step-size, in general.

The Malitsky-Mishchenko method involves, in addition to the step-size , another sequence {} which is (formally) initialized by <sup>0</sup> = +∞. Then, the scheme proceeds by updating a running estimate of the inverse Lipschitz constant of the gradient as the step-size

$$\alpha\_n = \min\left(\sqrt{1 + \theta\_{n-1}}\alpha\_{n-1}, \frac{||w^n - w^{n-1}||}{2||\nabla f(w^n) - \nabla f(w^{n-1})||}\right) \tag{3.21}$$

and gradient descent steps (3.19), see Alg. 1 in Malitsky and Mishchenko (2019). The sequence is updated in each step using the step-size = */*−1. The method requires an initial step size 0, which we take as for the gradient descent scheme (3.20).

#### **3.3.3 The Barzilai-Borwein method**

Originally devised as a special type of Quasi-Newton method, the scheme proposed by Barzilai and Borwein (1988) may also be interpreted as an adaptive gradient descent method using the step size

$$\alpha\_n = \frac{||w^n - w^{n-1}||^2}{(\nabla f(w^n) - \nabla f(w^{n-1}))^T (w^n - w^{n-1})}.\tag{3.22}$$

For the Barzilai-Borwein method also take the step-size of the gradient descent scheme (3.20) as 0.

The Barzilai-Borwein (3.22) leads to -linear convergence when applied to strongly convex quadratic objective functions (Dai and Liao, 2002). By a perturbation argument (Liu and Dai, 2001), the method can also be shown to be locally -linearly convergent when applied to general strongly convex and twice differentiable objective functions. However, there are counter-examples to global convergence for the latter class, leading to the development of non-monotone backtracking techniques (Fletcher, 2005) and step-size limiting approaches (Burdakov et al., 2019), for instance. When applying the Barzilai-Borwein method to computational micromechanics (Schneider, 2019), these globalization techniques were not necessary for the computational examples considered. On the contrary, in most cases they decreased the performance profile. For computing Laguerre tessellations, we made a similar observation, i.e., we did not encounter a single case where the Barzilai-Borwein method diverged, see section 3.4 below. Thus, we shall rely upon the Barzilai-Borwein prescription (3.22) *without* globalization, although we are not aware of a mathematical justification for this choice.

#### **3.3.4 Newton's method**

For finding a critical point of a twice continuously differentiable function , Newton's method iteratively updates

$$w^{n+1} = w^n + \alpha\_n \xi^n,\tag{3.23}$$

where the Newton increment solves the equation

$$
\nabla^2 f(w^n) \,\xi^n = -\nabla f(w^n).\tag{3.24}
$$

 ∈ (0*,* 1] is a damping factor (similar to the step-size of the gradient method). If \* is a critical point of with non-degenerate Hessian ∇2( \* ) and ∇2 satisfies a Lipschitz condition close to \* , Newton's method with ≡ 1 can be shown to be quadratically convergent in the vicinity of \* (Kantorovich, 1948). However, global convergence is not achievable that way. Among the various globalization techniques for Newton methods (Nocedal and Wright, 1999), the damped Newton method (3.23) is particular popular. Far away from a critical point, *<* 1 may be necessary to ensure global convergence. However, in the vicinity of a critical point, no damping is necessary ( ≡ 1 for ≥ 0). In this way, the local convergence speed of Newton's method is preserved.

If is (strictly) convex, determined by (3.24) is always a descent direction for at . Thus, a proper backtracking strategy allows us to determine a suitable step size . More precisely, given a backtracking factor ∈ (0*,* 1), starting from 0 = 1, the step size is iteratively shrunk, i.e., +1 = , until actual descent is observed, either in function value

$$f(w^n + \alpha\_n^k \xi^n) < f(w^n)$$

or in its gradient

$$\left\| \nabla f(w^n + \alpha\_n^k \xi^n) \right\| < \left\| \nabla f(w^n) \right\|.$$

For the problem at hand (3.17), we shall rely upon the latter criterion, as ∇ is cheaper to compute than . Under suitable conditions, this backtracking strategy can be shown to be globally convergent under technical conditions, see Boyd and Vandenberghe (2004). We take = 1*/*2 as our backtracking factor.

In the context of Laguerre tessellations with prescribed volumes, a further issue needs to be addressed. For computing the Newton increment , the equation (3.24) needs to be solvable. In section 3.2.1, we pointed out that constant vectors are in the kernel of ∇2() for any ∈ R . However, this is no problem if equation (3.24) is interpreted as an equation on , i.e., ∇( ) ∈ is given and ∈ is sought. A simple remedy consists in adding the (suitably scaled) projector onto the kernel onto ∇2( ), see equation (3.9).

Even restricted to , the function is not globally strictly convex. At such a point , also the Hessian will be degenerate (even when restricted to ). More precisely, taking a look at the Hessian of , we see that an empty Laguerre cell implies that both the -th row and the -th column only contain zeros.

A workaround has been proposed by Lévy (2015) by adding the constraint of non-empty cells to the backtracking strategy of Newton's scheme. This approach has been analyzed in Kitagawa et al. (2019) and shown to be globally convergent when started at a feasible point (i.e., with non-empty Laguerre cells).

However, for grain-size distributions with large variation, we found this cautious backtracking strategy to be computationally inefficient. Although the method is quadratically convergent locally, the first few steps may require extensive backtracking. As a simple alternative, we add - in case of vanishing -th cell - a suitably scaled 1 on the -th diagonal

entry of the Hessian. In this way, we ensure obtaining a descent direction and avoid the extensive backtracking.

#### **3.3.5 BFGS**

For large number of grains, the computationally most demanding step of Newton's method (3.23) is solving the linear system (3.24) which determines the Newton increment. If a direct solver is used, on the order of <sup>3</sup> floating point operations are required, which may be excessive for large .

To reduce the effort, Quasi-Newton methods approximate the inverse of the Hessian ∇() <sup>−</sup><sup>1</sup> directly, see Nocedal and Wright (1999). The most popular Quasi-Newton method is the Broyden-Fletcher-Goldfarb-Shanno method (Broyden, 1970; Fletcher, 1970; Goldfarb, 1970; Shanno, 1970), or BFGS, in short. In a nutshell, for every , the next iterate is computed as for Newton's method (3.23)

$$w^{n+1} = w^n + \alpha\_n \xi^n,$$

but the increment is determined by

$$
\xi^n = -H\_n \nabla f(w^n),
$$

where ∈ R × is a sequence of approximations to the inverse Hessian, which is, for BFGS, updated via

$$H\_{n+1} = H\_n + \frac{s\_n^T y\_n + y\_n^T H\_n y\_n}{(s\_n^T y\_n)^2} s\_n s\_n^T - \frac{1}{s\_n^T y\_n} \left( H\_n y\_n s\_n^T + s\_n y\_n^T H\_n \right) \tag{3.25}$$

where = ∇( +1) − ∇( ) and = +1 − . For strictly convex objective function , the BFGS update (3.25) preserves symmetry and positive definiteness of the Hessian approximation, see Nocedal and Wright (1999). Thus, if the initial approximation <sup>0</sup> is symmetric positive definite, all further will remain symmetric and positive definite. The effectiveness of BFGS crucially depends on the initial choice 0. A common choice (Nocedal and Wright, 1999) is to use the Barzilai-Borwein step size (3.22), i.e.,

$$H\_0 = \frac{\|w^0 - w^{-1}\|^2}{(\nabla f(w^0) - \nabla f(w^{-1}))^T (w^0 - w^{-1})} \operatorname{Id},$$

provided some inital gradient step has been taken, as for the Barzilai-Borwein method. When enriched by a backtracking strategy for determining the step size in (3.23), BFGS can be shown to converge for strongly convex functions (and further technical conditions) with a superlinear rate (Nocedal and Wright, 1999).

Notice that we do not need to pay specific attention to the kernel of ∇2 consisting of component-wise equal weights. This is a result of approximating the inverse of ∇2 on . Indeed, <sup>0</sup> is a multiple of the identity (with a proper scaling) and the BFGS update (3.25) only modifies .

The BFGS method reduces the complexity required for the Newton step from (<sup>3</sup> ) to (<sup>2</sup> ) (due to the update of and the matrix-vector product required for computing ). To further reduce the complexity to (), limited-memory BFGS ( − ), see Nocedal (1980), may be used. In general, however, the superlinear convergence is lost. Also, −  is often inferior to the Barzilai-Borwein method, as the latter scheme's non-monotonicity can be advantageous (Schneider, 2019; Wicht et al., 2020b).

## **3.4 Numerical demonstrations**

### **3.4.1 Setup and used hardware**

The algorithms described in sections 3.2.2 and 3.3 were implemented in Python 3.7 with Cython extensions. For computing Laguerre tessellations, we rely upon the Voro++ (Rycroft, 2009) code which we exposed to Python via Cython bindings using an OpenMP parallelization.

To reduce the stochastic influence, we use the Sobol low-discrepancy sequence (Sobol, 1967) for drawing log-normal grain-size distributions (Johnson et al., 1994). The numerical experiments were conducted on desktop computer with 6 3*.*7 GHz CPUs and 32 GB RAM.

For solving the linear system for the Newton increment (3.24), we rely upon NumPy's linalg.solve routine, which wraps LAPACK (Anderson et al., 1999). The "higher-level" vector and matrix operations are realized by NumPy, for instance the BFGS update (3.25).

For all convex optimization methods of section 3.3, we initialize the weights to zero. Using the estimate of Lyckegaard et al. Lyckegaard et al. (2011), as suggested in Bourne et al. Bourne et al. (2020), did not decrease the run-time, on average.

To ensure comparability between different solvers, all experiments are initialized with identical seeds. We realize this by *using* the 3D Sobol sequence on the (scaled) unit cube. The first few elements of the 3D Sobol sequence are a little pathological, as they correspond to the corners of the cell and some other very regularly positioned spots, leading to *extremely regular* Laguerre tessellations. Therefore, we skip a few elements of the Sobol sequence. To ensure comparability, we fix the number of elements to be skipped as 1234 (for no specific reason).

At this point, we wish to draw attention to the convergence criteria and tolerances used. For the convex optimization solvers, we test convergence by

$$\frac{\|\nabla f(w)\|}{\min\_{i=1}^{K} \phi\_i} \le \text{tol} \tag{3.26}$$

with tol = 10<sup>−</sup><sup>2</sup> . This choice ensures that the volume fraction of the smallest grain is accurate to two significant digits (which we assume sufficient for materials-science purposes). For centering, we use equation (3.11), i.e.,

$$\sqrt{\frac{\sum\_{i=1}^{K} \phi\_i \, d(x\_i^k, c\_i^k)^2}{(\max L\_i)^2}} \le \text{tol}$$

with tol is 10<sup>−</sup><sup>1</sup> . As Bourne et al. (2020), we noticed no visible improvement of the tessellation quality for lower values of tol. In contrast, the increase in computation times became quite visible.

We abbreviate grain-size distribution as GSD throughout section 3.4. For all the experiments in this section, we use cubical cells, i.e., <sup>1</sup> = <sup>2</sup> = 3. The prescribed algorithms work for anisotropic unit cells ("bricks") as well.

### **3.4.2 Computing the weights of a Laguerre tessellation with prescribed volume fractions**

The overall goal of this section is to evaluate the solvers introduced in section 3.3. For that purpose, we evaluate their iteration counts and run-times for problems of different complexity and grain count.

As our first example, we seek Laguerre tessellations with a mono-sized GSD, i.e., for grains, each grain should have a volume fraction 1*/*. As our smallest number of grains, we choose = 250. For smaller grain count, the run-times become even more negligible and a sensible evaluation is difficult. Starting from = 250, we subsequently double the grain count up to 16000 cells.

Larger numbers are possible (and computationally feasible), but probably limited in terms of applicability to computational homogenization (Shantraj et al., 2015; Wicht et al., 2020a).

The iteration counts of the investigated solvers are listed in Tab. 3.1. The gradient method (3.19) serves as our basic reference. Even for = 250, more than 300 iterations are required. For higher , the iteration counts are proportional to , i.e., they roughly double when doubling the number of grains.

For the Malitsky-Mishchenko method (3.21), the iteration counts are significantly lower than for the gradient method, roughly by one or two orders of magnitude. Still, we see a dependence of the iteration count on the grain number. More precisely, the iteration count depends approximately linearly on , but with a non-trivial offset at = 0.

The Barzilai-Borwein method leads to a consistently lower iteration count than the Malitsky-Mishchenko method, except for = 250, but no more than a factor of 2. BFGS (3.25) outperforms the Barzilai-Borwein scheme consistently in terms of iteration counts, also by no more than a factor of two.

Newton's method, see section 3.3.4, operates on a completely different level in terms of iteration counts, requiring 2 − 3 iterations for all grain counts considered. In particular, we see quadratic convergence in action.


**Table 3.1:** Iteration counts for different solvers for mono-sized GSD.

When evaluating computational methods, apart from the iteration counts, also the run-times are of particular importance, see Tab. 3.2. Even at first glance we see that the gradient method is not competitive in terms of run-time, which is certainly not surprising. The other four solvers perform quite well in terms of run-time for the majority of problems, requiring only a few seconds.

The Malitsky-Mishchenko and the Barzilai-Borwein method are both gradient methods of similar type, i.e., the computational expense of their individual steps are comparative (see also Tab. 3.3). For ≤ 2000, both of their run-times are similar (and similarly negligible). Only starting at = 4000, differences are visible. For = 16000, the Barzilai-Borwein method is about twice faster than Malitsky-Mishchenko.

Jumping ahead to Newton's scheme for the time being, we see that for small (up to = 2000), Newton's method is the fastest method investigated. For larger number of grains, Newton's scheme falls back behind the Barzilai-Borwein method in terms of run-time.

BFGS appears to be in between Barzilai-Borwein and Newton. For small grain count, it is slower than Newton's method, and for large , it is slower than the Barzilai-Borwein scheme.


**Table 3.2:** run-time in s for different solvers for mono-sized GSD.

More insight into these observations may be gained by investigating the run-time per iteration of the different solvers under consideration, see Tab. 3.3. Remarkably, the three gradient-type methods require

roughly the same computational effort per iteration. The expense is slightly higher for the Barzilai-Borwein method than for the other two gradient-type methods because the Python overhead plays a more distinct role, essentially due to the lower number of iterations required for Barzilai-Borwein method. Recall that we exposed the C++ - library Voro++ (Rycroft, 2009) to Python via Cython, and the major computational workload of the gradient-type methods amounts to computing Laguerre tessellations via Voro++. This library computes the Laguerre cells individually (and sequentially), using cell-linked lists (Allen and Tildesley, 1987). As our initial seeds are more or less equally distributed, and - for mono-sized GSD - we look for cells of equal volume, the effort for computing a single Laguerre cell is effectively independent of the cell count. In turn, the total computational expense for computing the entire Laguerre tessellation is proportional to . We see this trend quite clearly in Tab. 3.3 for the gradient-type solvers, highlighting the quality of the Voro++ - implementation. Notice the non-perfect scaling of the run-times, which is a result of the Python overhead and becomes less pronounced for large . A close look at the run-time per iteration for BFGS reveals a quadratic dependence on . In contrast, the effort required for a Newton step grows cubically in . These observations also clarify why the overall run-times of BFGS lag behind those of Newton's method for large . BFGS has a lower effort per iteration than Newton's method, but the overall iteration count is also significantly larger than for Newton's method, tipping the overall scale towards Newton.


**Table 3.3:** run-time per iteration in 10−<sup>2</sup> s for different solvers and a mono-sized GSD.

To finish this first example, we closer examined the convergence behavior of the solvers under consideration. For that purpose, for = 1000 grains, we solve equation (3.2) up to a precision of tol = 10<sup>−</sup><sup>10</sup>. This level of accuracy is certainly beyond what is necessary for applications, but gives interesting insights into the internal mechanisms of the solvers. The residual (3.26), is plotted vs. the iteration count in Fig. 3.2. We did not include Newton's method due to its quadratic convergence. It appears that all shown solvers converge linearly. The gradient method converges slowly, but in a steady and monotonic fashion. Both the Malitsky-Mishchenko (MM) and the Barzilai-Borwein (BB) method do not exhibit monotone convergence behavior but lead to an oscillation of the residual by one or two orders of magnitude. Overall, both schemes converge -linearly. This complies with theoretical estimates, see section 3.3. Overall, the convergence speed of BB is higher than for MM. BFGS also converges linearly, but in a monotonic fashion. Put differently, we have -linear convergence. BFGS exhibits a higher convergence rate than the other solvers shown. However, no superlinear convergence is observed. This is consistent to theory, as to enable superlinear convergence, the approximation to the inverse Hessian needs to be close, as well, see section 8.4 in Nocedal and Wright (1999). As the Hessian is a 1000 × 1000-matrix, the less than 100 iterations (and in turn, the corresponding matrix updates (3.25)) are simply insufficient to match the true inverse Hessian closely.

**Figure 3.2:** Convergence behavior of different inner solvers for the mono-sized GSD and = 250 grains.

As our next computational setup we consider a grain-size distribution with exactly two (distinct) volume fractions, each of which is equally probable. This may be viewed as a simplified bimodal grain-size distribution (Flipon et al., 2020a). Following Bourne et al. (2020), we use a factor of 5 between the two respective volume fractions. Thus, we consider grains with volume fraction 1*/*3 and grains with volume fraction 5*/*3, each. To keep consistency to the mono-sized GSD, the total number of grains 2 is varied between 250 and 16000, incrementally by factors of 2. Due to its lack of competitiveness for the mono-sized GSD, we did not consider the gradient descent scheme any further.


**Table 3.4:** Iteration counts for different solvers and the two-sized GSD.

The iteration counts for the four investigated solution methods are listed in Tab. 3.4. Both for Malitsky-Mishchenko and for Barzilai-Borwein, the iteration count increases monotonically with . A similar trend may be observed for Newton's method and BFGS. However, there are exceptions between = 250 and = 1000. It is also evident that the iterations counts for the two-sized GSD are consistently higher than for the mono-sized GSD, see Tab. 3.1.

Comparing the iteration counts of Malitsky-Mishchenko and Barzilai-Borwein, we see that for increasing , the discrepancy between the two solvers, in terms of iteration count, is increasing. For = 8000, Malitsky-Mishchenko requires about three times the iteration count of Barzilai-Borwein.

As for the mono-sized GSD, BFGS requires less iterations than Barzilai-Borwein, roughly by a factor of two for the larger .

For Newton's method, no quadratic convergence behavior can be observed for the two-sized GSD. Indeed, up to 7 Newton iterations are necessary to trigger the convergence criterion. Intuitively, this appears clear. Newton's method exhibits its strength close to the minimum. However, the seeds are evenly distributed, and the weights are equal initially. Thus, initially, the computed Laguerre cells will have roughly equal volume. In particular, they do not appear to serve as a promising starting value for Newton's method. Therefore, a few preliminary descent steps are necessary to bring the iterates into the domain of quadratic convergence.


**Table 3.5:** run-time in s for different solvers and the two-sized GSD.

The required run-times are listed in Tab. 3.5. In general, the run-times exceed those of the mono-sized problem, see Tab. 3.2. The Malitsky-Mishchenko method consistently lags behind the Barzilai-Borwein solver in terms of run-time. Newton's method is very strong except for the highest grain size considered. More precisely, up to = 2000, Newton's method outperforms all other solvers considered. For = 4000 and = 8000, the Barzilai-Borwein method is fastest. Again for the twosized scenario, BFGS does not perform admirably.

Notice that up to = 2000, the run-times of Barzilai-Borwein and Newton's method are minimal (up to about 2s). Also, for the highest grain count considered, 16000 grains, the Barzilai-Borwein solver needs less than 15s. Compared to the mono-sized problem, see Tab. 3.2, this amounts to an overall increase by a factor less than two.

For further analysis, we discard BFGS.

To model more realistic grain-size distributions, we use a log-normal distribution for the equivalent radius of the grains. To be more precise, for a grain of volume , we associate the radius of a sphere with precisely this volume to it

$$V = \frac{4\pi}{3}r^3, \quad \text{ie}, \quad r = \sqrt[3]{\frac{3V}{4\pi}}.$$

It is empirically known (Spettl et al., 2014) that the distribution of this equivalent radius follows a log-normal distribution with probability density function (Johnson et al., 1994)

$$\rho(r) = \frac{1}{r} \cdot \frac{1}{\sigma \sqrt{2\pi}} \exp\left(-\frac{(\ln\left(r\right) - \mu)^2}{2\sigma^2}\right),\tag{3.27}$$

depending on the two real parameters and . The mean and standard deviation compute as (Johnson et al., 1994)

$$\text{mean} = e^{\mu + \frac{1}{2}\sigma^2} \quad \text{and} \quad \text{stdev} = e^{\mu + \frac{1}{2}\sigma^2} \sqrt{e^{\sigma^2} - 1}.$$

As the mean equivalent radius plays no role in the geometry generation, we set it equal to unity. In turn, we obtain the relation = −1*/*2 2 , reducing the original two parameters to the single parameter . In turn, the expression for the standard deviation simplifies to

$$\text{stdev} = \sqrt{e^{\sigma^2} - 1}.\tag{3.28}$$

Thus, for given standard deviation stdev, we may determine the parameter characterizing the log-normal distribution (3.27) with mean 1 explicitly

$$
\sigma = \sqrt{\log(1 + \text{stdev}^2)}.
$$

Thus, we may parameterize grain-size distributions by the standard deviation. For vanishing standard deviation, we recover the mono-sized distribution. Spettl et al. (2014) report stdev = 0*.*35 as typical for a polycrystalline material that underwent grain-boundary migration by capillary effects.

To sum up, we prescribe the unit cell size (1*,* 2*,* 3), the number of grains and the standard deviation stdev of the log-normal distribution. Then, we sample radii 1*,* 2*, . . . ,*  according to the log-normal distribution by use of the Sobol sequence (Sobol, 1967). Then, we compute the appropriate volume fractions by

$$\phi\_i = \frac{r\_i^3}{\sum\_{j=1}^K r\_j^3},$$

which serves as input for our optimization routines.

For the experiment at hand, we fix the grain count to = 1000 and vary the standard deviation from 0 to 0*.*5 in 0*.*05-steps. Fig. 3.3 contains images for microstructures generated for varying standard deviation. The seeds are identical for all structures. We see that some grains grow, whereas others shrink, for increasing standard deviation. Please note that the colors appearing in the figures of this article showing

polycrystalline microstructures are chosen randomly. Their sole purpose is to distinguish the individual cells. In particular, they are devoid of any physical meaning.

**Figure 3.3:** Influence of increasing standard deviation for polycrystalline microstructures with log-normally distributed equivalent radii and 1000 grains. The seed locations are fixed and the plots show different draws from the log-normal distribution. (a) stdev = 0*.*00, (b) stdev = 0*.*10, (c) stdev = 0*.*20, (d) stdev = 0*.*30, (e) stdev = 0*.*40, (f) stdev = 0*.*50.

The iteration counts and run-times for increasing standard deviation and the Malitsky-Mishchenko and the Barzilai-Borwein scheme as well as Newton's method are listed in Tab. 3.6.


3 Fast methods for computing Laguerre tessellations with prescribed volume fractions

**Table 3.6:** Iteration counts and run-times for different solvers and varying standard deviation for = 1000 cells.

Apparently, as a general rule of thumb, increasing standard deviation also leads to an increase in iteration count *and* run-time for all solvers considered. This is not only caused by the increased difficulty, but also by the convergence criterion (3.26) in use. Indeed, the smallest among the volume fractions enters the convergence criterion. As the smallest volume fraction decreases for increasing standard deviation, the convergence criterion becomes stricter, as well. At this point, we may question the convergence criterion we use. However, this convergence criterion is standard, see, for instance, Lévy (2015).

Taking a closer look at the iteration counts in Tab. 3.6, we see that the Malitsky-Mishchenko method consistently requires more iterations than the Barzilai-Borwein method. For stdev = 0*.*5, it even requires more than thrice the iteration count of Barzilai-Borwein. The required iterations for the Barzilai-Borwein method increases superlinearly in the standard deviation. It appears that the iteration count roughly doubles every 0*.*1-step in standard deviation.

For increasing standard deviation, Newton's method is still superior, but with a smaller margin, because it can profit less from its local quadratic convergence behavior. For low standard deviation (below 0*.*15), three iterations suffice. For intermediate standard deviation (less than 0*.*4), only a few Newton steps need to be taken before entering the domain of quadratic convergence. For larger standard deviation, the iteration counts increase significantly, rising up to 48 for stdev = 0*.*5.

Recall that we used = 1000 grains. For stdev = 0, we recover the mono-sized experiment. The run-times for this case lie well below one second. For Malitsky-Mishchenko, compared to the mono-sized case, for stdev = 0*.*5 the run-time is increased by a factor of 100. For the other two solvers, the increase is less severe, but still noticeable. Overall, for standard deviations up to 0*.*4, Newton's method is extremely fast, with run-times below one second. Only for higher standard deviation, a significant increase in run-time is observed. The run-times for the Barzilai-Borwein method increase more steadily with increasing standard deviation, and are only slightly above those of Newton's method. Recall, from the previous experiments, that the optimal solver strongly depends on the number of grains present in the tessellation. For = 1000 grains, Newton's method always proved the fastest. For the last experiment, we saw that the standard deviation of the grain-size distribution also plays a significant role when evaluating the performance of the solvers in question.

#### **3.4.3 Computing centroidal Laguerre tessellations for prescribed volume fractions**

In the previous section 3.4.2, we computed periodic Laguerre tessellations for prescribed volume fractions and *fixed* seeds. In this section, we wish to determine centroidal Laguerre tessellations for prescribed volume fractions, and add the centering algorithms of section 3.2.2 on top of the optimization algorithms of section 3.3. We have seen that either Newton's method or the Barzilai Borwein scheme proved computationally most efficient for the problems in section 3.4.2. Thus, we restrict

to these two solvers for the subsequent investigations. For the entire section, we use the convergence criterion (3.11) with tolerance tol = 10<sup>−</sup><sup>1</sup> . Recall that we check convergence for the convex optimization solvers by condition (3.26) with tol = 10<sup>−</sup><sup>2</sup> .

Consistent to section 3.4.2, we first investigate a mono-sized GSD with increasing number of grains.

The number of centering iteration counts is listed in Tab. 3.7. We see that the number of centering iterations depends only on the "outer" algorithm used, and is independent on the "inner" algorithm. This fact appears clear at first, but becomes more interesting as we chose the tolerance for inner residual (3.26) comparatively large. Indeed, for inexact Newton methods (Dembo et al., 1982), the choice of the inner tolerance has a strong influence on the overall performance of the scheme, see Dembo et al. (1982). However, this seems not to be the case for the problem at hand.


**Table 3.7:** Centering iteration counts for different optimization solvers, combined with the Lloyd-type algorithm Alg. 1 and its Anderson-accelerated version Alg. 2, and the mono-sized GSD.

Comparing the Lloyd-type algorithm Alg. 1 and its Anderson-accelerated version Alg. 2, we see Anderson acceleration decreases the iteration count in any case. However, for high grain count, the improvement is not large. The associated run-times are listed in Tab. 3.8. Up to = 2000 grains, Newton's method is faster. For higher grain count, the Barzilai-Borwein method leads to a lower run-time. Comparing the centering algorithms, Anderson acceleration consistently lowers the run-times.


**Table 3.8:** Run-time in s for different optimization solvers, combined with the Lloyd-type algorithm Alg. 1 and its Anderson-accelerated version Alg. 2, and the mono-sized GSD.

More insight might be gained from Fig. 3.4, where the number of inner iterations is shown vs. the current centering iteration for = 250. For Lloyd's method, see Alg. 1, the inner iteration count decreases monotonically for increasing centering iteration. This property is preserved by Anderson acceleration (except for a single centering iteration).

To sum up, for large grain counts, i.e., ≥ 4000, combining the Barzilai-Borwein method with Alg. 2 performs best, requiring less than 20s for the largest problem considered.

**Figure 3.4:** Inner iterations vs centering iterations for the outer solvers for the mono-sized GSD and = 250 grains.

As our next example, as in section 3.4.2, we consider a two-sized grainsize distribution with a ratio of 5 between the volume fractions of the two "phases", each of which has grains.

The centering iteration counts are listed in Tab. 3.9. Compared to the mono-sized GSD, the iteration counts are slightly higher, on average. Still, the observations for the mono-sized GSD also apply for this twosized case.


**Table 3.9:** Iteration counts for different optimization solvers, combined with the Lloyd-type algorithm Alg. 1 and its Anderson-accelerated version Alg. 2, and the two-sized GSD.

As already observed in Bourne et al. (2020), the time per centering iteration decreases. Therefore, we investigate the corresponding runtimes, contained in Tab. 3.10. Already for = 4000, Barzilai-Borwein combined with Alg. 2 outperforms all other variants. For the highest grain count considered, the latter combination converges in about half a minute, which is about 33% higher than for the mono-sized GSD.


**Table 3.10:** run-time in s for different optimization solvers, combined with the Lloyd-type algorithm Alg. 1 and its Anderson-accelerated version Alg. 2, and the two-sized GSD.

As our last example in this paragraph, we consider the problem of computing a centroidal periodic Laguerre tessellation with = 1000 grains with log-normally distributed GSD, cf equation (3.27), and varying standard deviation.

The number of centering iterations, see Tab. 3.11, is almost independent of the inner solver (different only by 1, at most). We see that the number of centering iterations increases monotonically with increasing standard deviation (except for Anderson and the transition from 0*.*45 to 0*.*5). However, in contrast to the inner iterations, see Tab. 3.6, the increase is modest. For both centering algorithms, the run-times are comparable, and thus mainly determined by the inner solver. Taking a closer look at the run-times for Lloyd-type centering, we see that Newton's method serves as the faster inner solver up to a standard deviation of 0*.*2. However, Barzilai-Borwein is only slightly faster. Indeed, for stdev = 0*.*5, both inner solvers lead to an overall run-time less than half a minute.


3 Fast methods for computing Laguerre tessellations with prescribed volume fractions

**Table 3.11:** Centering iteration counts and run-times for different optimization solvers, combined with the Lloyd-type algorithm Alg. 1 and its Anderson-accelerated version Alg. 2, and varying standard deviation. \* means that the method did not converge in a reasonable amount of time.

The generated centroidal polycrystalline microstructures, for different values of the standard deviation, are shown in Fig. 3.5. For comparison, the non-centroidal versions may be of interest, see Fig. 3.3. Recall that we used identical seeds for all cases, i.e., independent of standard deviation and whether centering is applied or not. Also, the coloring schemes are kept consistent. Thus, it is possible to match the various grains for all cases considered. Apparently, the centroidal tessellations in Fig. 3.5 exhibit a much higher shape regularity than the non-centered tessellations of Fig. 3.3. Also, the sphericity of the large grains becomes more pronounced for increasing standard deviation.

**(d) (e) (f)**

**Figure 3.5:** Influence of increasing standard deviation for *centered* polycrystalline microstructures with log-normally distributed equivalent radii and 1000 grains. (a) stdev = 0*.*00, (b) stdev = 0*.*10, (c) stdev = 0*.*20, (d) stdev = 0*.*30, (e) stdev = 0*.*40, (f) stdev = 0*.*50.

## **3.5 Examples with a higher level of complexity**

After thoroughly assessing the performance of the convex-problem solvers and the centering, we shall investigate examples with a higher degree of complexity.

We used the Barzilai-Borwein method to generate a microstructure with 100000 = 10<sup>6</sup> grains of identical volume, up to an accuracy of 10<sup>−</sup><sup>2</sup> , as before. For quasirandom initial positions, the Barzilai-Borwein method needed 282*.*85s, i.e., less than 5 minutes. The resulting microstructure is shown in Fig. 3.6a. Coupled to the Anderson-accelerated Lloyd algorithm, solved up to an accuracy of 10<sup>−</sup><sup>4</sup> , took 338*.*32s, i.e., less than six minutes. The resulting microstructure is shown in Fig. 3.6b and may be compared to its non-centered cousin, cf Fig. 3.6a. Even without the multi-level strategies of Mérigot (2011) and Lévy (2015), we are able to generate microstructures with high complexity to be used in crystal-plasticity FEM (Becker, 1991; Barbe et al., 2001; Roters et al., 2010) or FFT (Lebensohn and Rollett, 2020; Shantraj et al., 2015; Wicht et al., 2020a) frameworks in a matter of seconds to a few minutes.

**Figure 3.6:** Generated microstructure with 100000 grains of identical volume. (a) Noncentroidal Laguerre tessellation, (b) Centroidal Laguerre tessellation.

The framework presented may also be used for generating microstructures with a morphological texture, i.e., a prescribed spatial gradient in the grain sizes. Such microstructures may be the result of high temperature gradients that may occur inside of polycrystalline metals, for example during welding (Lehto et al., 2016).

We follow the ideas of Bourne et al. (2020) and restrict the initial positions of the seeds to particular subsets. Of course, the resulting grains need not be entirely contained in the pre-selected areas, in particular if coupled to the centering methods. Still, interesting microstructure emerge, as we shall demonstrate by an explicit example.

We consider a cubical cell [0*,* ] × [0*,* ] × [0*,* ]. We shall divide the grains into three classes:


Thus, we have three types of grains. Each type has the same total volume. However, the grain volume of class one is four times smaller than for the second class, and sixteen times smaller than the third class. Also, the first class is restricted to the first quarter of the -face, whereas the second class is restricted to the second quarter, and the third phase makes up the last half of the -axes.

Generating the microstructure with Barzilai-Borwein and Anderson-Lloyd (with accuracy as above) took 2*.*0s. The resulting microstructure is shown in Fig. 3.7b. The first class is shown in red-to-black, the second class in yellow, whereas the third class is characterized by a grayish coloring.

**Figure 3.7:** Generated polycrystalline microstructure with spatial size-gradient. (a) Slice through the center, (b) Volumetric view.

As discussed in the introduction of this article, modern image-based characterization techniques enable measuring grains-size distributions of polycrystalline materials empirically. Although analytically given grains-size distributions have their merits, relying upon empirically given GSDs permits a close comparison to real-world specimens, in particular concerning their inherent imperfections and deviations from theoretical GSDs.

To demonstrate the capabilities of the proposed methodology, we use the emprical GSDs determined by Döbrich et al. (2004) as our target GSDs for generating synthetic microstructures. More precisely, we use histogram bins for describing the desired GSD, see Fig. 3.8.

**Figure 3.8:** Comparison of targeted and reached bin frequencies for the GSD of Döbrich et al. (2004). (a) 250 grains, (b) 500 grains, (c) 1000 grains.

Different centroidal periodic Laguerre tessellations were generated with an increasing number of grains using Barzilai-Borwein in combination with the Anderson-Lloyd method. The resulting microstructures are shown in Fig. 3.9, and the run-times are collected in Tab. 3.12. We see that even the largest structure with 8000 grains required only slightly more than two minutes to compute.

**Figure 3.9:** Generated microstructures with empirical GSD of Döbrich et al. (2004). (a) 250 grains, (b) 500 grains, (c) 1000 grains, (d) 2000 grains, (e) 4000 grains, (f) 8000 grains.


**Table 3.12:** Run-times for different numbers of grains and the empirical GSDs of Döbrich et al. (2004) using BB and Alg. 2.

Of particular importance to us is the ability of the methodology to match any prescribed empirical GSD. For this reason, we take a look at Fig. 3.8, where both the target and the realized histograms are shown for 250, 500 and 1000 grains. For 250 and 500 grains, differences between the histograms are visible, although on a level probably within engineering accuracy. For 1000 grains, the target and the realized histograms are difficult to distinguish by eye.

To gain more insight, we plotted the absolute errors (in %) between the desired and targeted frequencies per bin, cf Fig. 3.10. We see a maximum error of about 0*.*5% for 250 grains, 0*.*3% for 500 grains and below 0*.*15% for 1000 grains. For a higher grain count, the error decreases further. More precisely, we see that the error decreases as 1*/*, where denotes the number of grains, i.e., the error reduces roughly by 50% when doubling the number of grains. This is no surprise, but a result of our quasirandom sampling of the grain sizes, which theoretically predicts a 1*/* decrease of error upon sampling (Asmussen and Glynn, 2007). In contrast, a purely random sampling only reduces the expected error by a factor of 1*/* √ . To gain some intuition into this difference, essentially to arrive at the same error for random sampling that we obtained with quasirandom sampling and 1000 grains would require 1000<sup>2</sup> , i.e., one million, grains, on average. Last but not least let us remark that the maximum error decreases linearly, but the bin-wise error might not. Indeed, taking a look at the second-smallest bin, the error only decreases for 2000 grains compared to 250 grains. However, the corresponding error is already small for 250 grains.

**Figure 3.10:** Absolute frequency error of the microstructure realizations of Fig. 3.9, as a result of the quasi-random sampling, compared to the target values of Döbrich et al. (2004). (a) 250 grains, (b) 500 grains, (c) 1000 grains, (d) 2000 grains, (e) 4000 grains, (f) 8000 grains.

## **3.6 Conclusion**

This work has been devoted to generating geometry-based microstructure models for polycrystalline materials in terms of Laguerre tessellations. According to Ball and Carstensen (2015) Laguerre tessellations describe polycrystalline microstructures with convex grains. Put differently, in order to describe a geometrically more complex polycrystalline microstructure, non-convex grains have to be considered.

In a recent work, Bourne et al. (2020) have shown that, for any given set of seed points and prescribed volume fractions, a corresponding set of Laguerre weights realizing a Laguerre tessellation with the prescribed volume fractions may be determined by solving a convex program. As a consequence, using block-box optimization methods, Bourne et al. (2020) were able to generate polycrystalline microstructures of industrial complexity in the order of minutes.

In this article, we retained the basic algorithmic framework of Bourne et al. (2020) for computing the weights of a periodic Laguerre tessellation for prescribed seed locations and volume fractions, and exploited the benefits of modern non-monotone methods of gradient type, i.e., the Malitsky-Mishchenko method and the Barzilai-Borwein scheme, in the context of semi-discrete optimal transport.

In the numerical experiments we saw that, for small to moderate grain count and low standard deviation of the grain-size distribution, Newton's method outperforms all other methods considered. For larger number of grains or for larger spread of the grain sizes, the Barzilai-Borwein method proved superior in terms of run-time. As the run-time for small- to medium-sized problems is small anyway, the Barzilai-Borwein method may serve as a general-purpose method for computing the weights of Laguerre tessellations.

To increase the physical realism, we also contributed to computing centroidal Laguerre diagrams of prescribed volume fractions by suggesting to combine the classical Lloyd's algorithm to Anderson acceleration, a

general-purpose acceleration scheme for fixed-point methods. For all problems considered, the Anderson-accelerated Lloyd's method proved superior to the classical Lloyd's method. Anderson acceleration might also prove more powerful if higher accuracy for centering is required, in particular when compared to BFGS-type centering (Liu et al., 2016; Hateley et al., 2015). However, in the context of polycrystalline microstructures, such a high level of accuracy appears not to be necessary. All in all, combining the improved convex optimization solver, the centering and an OpenMP-parallel implementation permits generating polycrystalline microstructures of industrial complexity and relevance in a matter of seconds to minutes, at most, reducing the run-time and increasing the number of grains compared to the work of Bourne et al. (2020). In particular, for grain-size distribution with large standard deviation of the grain sizes, it appears imperative to accept iterates with vanishing Laguerre cells. Classically, backtracking methods are used to avoid empty Laguerre cells during the iterative process. The proposed gradient methods avoid such backtracking - in fact, their performance decreases when backtracking is used. Furthermore, we also proposed a simple regularization strategy for Newton's method in case of vanishing Laguerre cell volume which leads to a robust and competitive computational scheme.

Our solution schemes may be further improved by using multi-level strategies, as pioneered by Mérigot (2011) and Lévy (2015). However, for applications to micromechanics, the number of grains to consider typically does not exceed a few thousand.

Notice that we describe the morphological texture of a polycrystalline aggregate in terms of the seed locations and the volume fractions, taking into account only one-point statistics (Torquato, 2002). We have shown, by example, that microstructures with a spatial gradient in the grain-size distribution may be generated. Also, experimentally determined grainsize distributions may be reproduced with high accuracy, provided the equivalent radii are drawn in a pseudo-random fashion. Still, it might

be interesting to see whether the algorithms could be extended to more general morphological features, e.g., sphericity.

Last but not least it is clear that, in order to conduct computational experiments, the generated polycrystalline microstructures need to be furnished with crystal orientations per grain, leading to possible crystallographic texture of the material under consideration. However, the latter is typically realized by a post-processing procedure, see, e.g., Quey et al. (2018), and is thus complementary to the approach presented.

### **Chapter 4**

# **Generating polycrystalline microstructures with prescribed texture coefficients <sup>1</sup>**

## **4.1 Introduction**

The methods proposed in Chap. 3 enable the generation of polycrystalline microstructures of industrial complexity within several seconds up to a few minutes. As outlined in Sec. 2.4 the next step in the generation of synthetic polycrystalline microstructure is assigning crystallographic orientations to each grain.

In this chapter, we propose a method based on condensing the information carried by the crystallite orientation distribution function (CODF) via the coefficients of a tensorial series expansion (Guidi et al., 1970; Adams et al., 1992). As these tensorial coefficients are easy to compute from experimental data, directly linked to the bounds of mechanical properties (Böhlke and Bertram, 2003; Böhlke and Lobos, 2014) and easily applied in microstructural sensitive design (Fullwood et al., 2010; Lobos and Böhlke, 2015), we prefer them to an approach involving spherical harmonics (Roe, 1965; Bunge, 1982). Moreover, Böhlke (2005) introduced

<sup>1</sup> This chapter is based on the publication by Kuhn et al. (2022) whereas minor typographical and formatting changes have been made for cohesion of the manuscript and the readers convenience

a method for estimating the crystallite orientation distribution function for a finite number of given texture coefficients and Junk et al. (2012) analyzed this approach in the context of maximum entropy moment problems. Böhlke (2006) derived a hierarchy of evolution equations for the texture coefficients under the Taylor-Voigt assumption of vanishing strain fluctuations on the microscale. Motivated by these observations, we propose to use a limited set of tensorial texture coefficients as the input of a microstructure generator, with the goal to obtain crystallographic orientations for each crystallite matching the prescribed texture coefficients. We follow a two-step procedure. First, the orientations are sampled randomly. In a second step, these orientations are corrected in terms of a gradient-based optimization technique, which ensures that the resulting texture coefficients match the prescribed ones up to a given tolerance. To introduce the *T*exture coefficient *O*ptimization for *P*rescribing orientations (TOP) method, we briefly revisit the notation of rotations in Sec. 4.2.1. In Sec. 4.2.2, we outline the problem formulation and the solution scheme. We investigate the capabilities of the proposed method, by comparing it to state-of-the-art algorithms in Sec. 4.3.

## **4.2 Background on modeling and optimization**

### **4.2.1 Representing the texture**

Describing the deformation behavior of a polycrystal, i.e., an agglomerate of crystals, by crystal plasticity requires taking the distinct orientation of each crystallite into account. For instance the stiffness tensor C in Hooke's law (2.31) depends on the crystal orientation. To describe the orientations in a polycrystalline material, for each crystallite we use a

proper orthogonal tensor

$$Q = \sum\_{i=1}^{3} g\_i \otimes e\_i,\tag{4.1}$$

where (1*,* 2*,* 3) and (1*,* 2*,* 3) represent the fixed orthonormal basis of the sample and the crystallite, respectively. Thus, the orientation of a crystallite is encoded by the rotation from the crystal coordinate system into the sample coordinate system. All orientation tensors are elements of the group of proper rotations in three dimensions, i.e., ∈ (3). In the following, we use the expressions rotation and orientation interchangeably.

For a given polycrystal, the orientation may be succinctly described in terms of the crystallite orientation distribution function (CODF) , a probability distribution on (3) w.r.t. the Haar Measure (normalized to unity). More precisely, the CODF is non-negative,

$$f(Q) \ge 0 \quad \forall \quad Q \in SO(3),\tag{4.2}$$

and normalized

$$\int\_{SO(3)} f(Q) \, dQ = 1. \tag{4.3}$$

For a (measurable) subset ⊆ (3) of orientations, the expression ∫︀ () computes the probability to find orientations contained in the set . Due to the invariance properties of the Haar measure (Gel'fand et al., 2018), the invariance property

$$\int\_{SO(3)} f(Q) \, dQ = \int\_{SO(3)} f(QQ\_0) \, dQ \quad \text{holds for all} \quad Q\_0 \in SO(3). \tag{4.4}$$

Moreover, the CODF reflects the underlying symmetries of the crystals forming the aggregate, i.e.,

$$f(Q) = f(QH^C) \quad \text{holds for all} \quad H^C \in S^C \subseteq SO(3), \tag{4.5}$$

where denotes the (discrete) symmetry group of the crystals, as the orientation states and correspond to the same physical orientation state of the crystallite. As a result of from a forming process, the sample itself may possess a certain symmetry, encoded by a symmetry group . This is reflected by the CODF in terms of the condition

$$f(Q) = f(H^S Q) \quad \forall \quad H^S \in S^S \subseteq SO(3),\tag{4.6}$$

where denotes the symmetry group of the sample. For the sake of readability, we will only consider the case of cubic crystals and triclinic sample symmetry in the following (Morawiec, 2003, Ch. 3). Our approach permits a straightforward extension to the general case with arbitrary crystal and sample symmetries, see Zheng and Zou (2001a;b). If the CODF is not uniform, i.e., () ̸= 1 for some ∈ (3), then the material will be said to possess a crystallographic texture.

Working with the full CODF is oftentimes impractical. Guidi et al. (1970) proposed a way to condense the encoded information. More precisely, any square integrable function may be expressed in terms of a tensorial Fourier series (Guidi et al., 1970; Adams et al., 1992)

$$f(Q) = 1 + \sum\_{i=1}^{\infty} \mathbb{V}'\_{\langle \alpha\_i \rangle} \cdot \mathbb{F}'\_{\langle \alpha\_i \rangle}(Q),\tag{4.7}$$

involving the tensorial Fourier, or texture, coefficients V ′ ⟨⟩ and the tensor functions F ′ ⟨⟩ (Böhlke, 2005; Fernández and Böhlke, 2019). The subscript denotes the number of linear independent harmonic cubic tensors of rank . Harmonic tensors, denoted by a prime (·) ′ , are completely symmetric and traceless, i.e., the relations

$$A'\_{ijkl} = A'\_{jikl} = A'\_{jilk} = \cdots , \quad \sum\_{i=1}^{3} A'\_{iikl} = 0 \tag{4.8}$$

reduce the number of degrees of freedom to 2 + 1 (Böhlke, 2005). The tensor functions F ′ ⟨⟩ () in equation (4.7) arise by rotating suitable reference tensors T ′ ⟨⟩ (Guidi et al., 1970)

$$\mathbb{F}'\_{\langle \alpha\_i \rangle}(Q) = Q \star \mathbb{T}'\_{\langle \alpha\_i \rangle},\tag{4.9}$$

where

$$Q \star \mathbb{T} = T\_{ij\dots p}(Qe\_i) \otimes (Qe\_j) \otimes \dots \otimes (Qe\_p) \tag{4.10}$$

denotes the Rayleigh product, i.e., the rotation of a tensor by . Without loss of generality and following Böhlke et al. (2010) and Dyck and Böhlke (2020), we normalize the Frobenius norm of the reference tensors to

$$\|\mathbb{T}'\_{\langle \alpha\_i \rangle}\| = 2\alpha + 1 \tag{4.11}$$

instead of ‖T ′ ⟨⟩ ‖ = 1 (as done by Guidi et al. (1970)). Interpreting the texture coefficients V ′ ⟨⟩ as a convex combination of normalized reference tensors of single crystal states (Fernández and Böhlke, 2019), we compute the tensorial texture coefficients in the case of discrete (experimental) orientation data as (Böhlke, 2006)

$$\mathbb{V}'\_{\langle \alpha\_i \rangle}(Q\_1, \dots, Q\_K) = \frac{1}{2\alpha\_i + 1} \sum\_{k=1}^K c\_k \, Q\_k \star \mathbb{T}'\_{\langle \alpha\_i \rangle},\tag{4.12}$$

where denotes the volume fraction of orientation among orientations. From equations (4.11) and (4.12) it follows that the Frobenius norm of all texture coefficients V ′ ⟨⟩ lies within the interval [0*,* 1]. Böhlke (2006) and Böhlke and Lobos (2014) use the norm to measure the degree

of anisotropy. For a completely uniform CODF, all texture coefficients vanish, whereas for the case of single crystals the norm of all texture coefficients is equal to one.

As we seek a compact representation of the CODF in the following, we will restrict to the texture coefficients up to rank six. As the cubic reference tensor of rank two is zero and, because of the cubic crystal symmetry, odd-rank reference tensors up to a rank of eight vanish (Guidi et al., 1970), we focus on the texture coefficients of rank four and six. For a precise overview and thorough discussion of texture coefficients in a more general context, the reader is referred to the work of Fernández and Böhlke (2019).

## **4.2.2 Texture coefficient optimization for prescribing orientations (TOP)**

To create digital representations of polycrystalline microstructures, it is not sufficient to solely match the grain morphology. In addition, we have to take the orientation state, i.e., the CODF, into account (Kocks et al., 2000). Many tools either rely on simple model CODFs to generate orientations (Quey et al., 2011; 2018), need (possibly vast) experimental data (Groeber and Jackson, 2014) or at least a representation of the complete CODF (Eisenlohr and Roters, 2008) to generate orientations for digital microstructures. In this section, we propose a method to generate orientations based on tensorial texture coefficients. For a given unit cell, subdivided into individual grains, our goal is to prescribe the orientation per grain in such a way that the resulting texture coefficients V ′ ⟨⟩ of the unit cell match the prescribed ones V¯′ ⟨⟩ up to a given tolerance tol, thus approximating the underlying CODF. To this end, we formulate our objective function as the difference in independent components between

the current and the desired texture coefficients

$$\ell(Q\_G) = \sqrt{\sum\_{\beta=4}^{\alpha\_{\text{max}}} \|\bar{\mathcal{V}}'\_{\langle\beta\rangle} - \mathcal{V}'\_{\langle\beta\rangle}(Q\_G)\|^2},\tag{4.13}$$

where = (1*, . . . ,* ) is a vector of orientations (one for each of the cells) and max denotes the maximum rank of the considered texture coefficients.

Starting from a randomly initialized set of orientations, we seek to minimize the objective function *ℓ*(). The objective function *ℓ* is continuously differentiable, and it natural to use gradient-based optimization techniques. Indeed, the objective function *ℓ* is actually a polynomial in the components of the individual grain orientations . In particular, the function *ℓ* is infinitely often differentiable.

The simplest conceivable approach proceeds via gradient descent, i.e., following the direction of steepest descent. Please keep in mind that the objective function *ℓ* is defined (4.13) on the -fold product of (3). The space of special orthogonal tensors (3) forms a non-linear subset of the space of second-order tensors. In particular, a simple explicit Euler discretization of gradient descent does not work, as the next trial point does not necessarily lead to vector whose components satisfy the constraints of (3) (Taylor and Kriegman, 1994; Mäkinen, 2008).

This is illustrated in Fig. 4.1a. On a curved space with Riemannian metric, the natural extension of straight lines are so-called *geodesics*, which emanate from a point in a specific (tangent) direction by parallel translation. On a (compact, matrix) Lie group with its natural Riemannian metric (the Killing form), following the geodesics may be computed in terms of the matrix exponential. In case of (3), this reduces to Rodrigues' formula

$$\exp\left(J(\omega)\right) = Id + \frac{\sin(\theta)}{\theta}J(\omega) + \frac{1-\cos(\theta)}{\theta^2}J^2(\omega),\tag{4.14}$$

describing a rotation around an axis by an angle and where we set =  as well as

$$J(\omega) = \begin{pmatrix} 0 & -\omega\_3 & \omega\_2 \\ \omega\_3 & 0 & -\omega\_1 \\ -\omega\_2 & \omega\_1 & 0 \end{pmatrix}. \tag{4.15}$$

The gradient descent scheme works as follows, where a step size  *>* 0 is fixed beforehand. Suppose the -th iterate = ( 1 *, . . . ,*  ) is given ( = 0*,* 1*, . . .*). Then, we investigate the function

$$\ell\_i(\omega) = \ell(Q\_1^i \exp\left(J(\omega\_1)\right), \dots, Q\_K^i \exp\left(J(\omega\_K)\right)), \tag{4.16}$$

where ≡ (1*, . . . ,* ) ∈ R <sup>3</sup>. Computing ∇*ℓ*(0) ∈ R <sup>3</sup> by the chain rule, gradient descent proceeds via

$$Q\_G^{i+1} = \left( Q\_1^i \exp\left( J(-t\left[\nabla\ell\_i(0)\right]\_1) \right), \dots, Q\_K^i \exp\left( J(-t\left[\nabla\ell\_i(0)\right]\_K) \right) \right) \tag{4.17}$$

**Figure 4.1:** (a) Gradient descent on a linear space vs. (b) descent along a geodesic (dashed line) on the manifold (3).

### **4.3 Computational investigations**

#### **4.3.1 Setup**

In the following sections, we wish to provide insights into the performance of the proposed TOP method. As a first step, we consider linear elastic material behavior, i.e., we investigate the effective stiffness. Secondly, we study cyclic stress-strain hysteresis.

To create the morphology of the microstructures under investigation, we use the algorithm described in Kuhn et al. (2020) to generate digital polycrystalline microstructures with prescribed volume fractions. For the morphology we consider two cases, a unique and log-normal grain size distribution (GSD). The former means that all grains have the same volume, i.e., = <sup>1</sup>*/*, where denotes the total number of grains in the volume element. Restricting to a unique grain size permits us to study the influence of the individual grain orientations exclusively. Due to their frequent occurrence in experiments (Spettl et al., 2014), we also investigate microstructures with an equivalent diameter following a log-normal grain size distribution with mean equal to unity and a standard deviation of 0*.*15, see Kuhn et al. (2020).

For a fixed morphology, we furnish the grains with orientations, where we investigate three different CODFs, namely a uniform, one with a slight texture and one with an increased texture. For the former, we compare the accuracy of TOP and the algorithm proposed by Quey et al. (2011; 2018), integrated into the polycrystal generation software Neper. In addition, we include a random orientation sampling (realized using the scipy implementation for sampling the Haar distribution (Stewart, 1980)) as a benchmark. For the textured CODFs, we compare the TOP method to random sampling from discrete orientation measurements, which is a common practice (Gillner and Münstermann, 2017; Tu et al., 2019). For the textured CODFs and a log-normal GSD, we additionally consider the Texture Discretization Technique (TDT) algorithm proposed

by Melchior and Delannay (2006), which, in a first step, samples orientations using the method proposed by Tóth and Van Houtte (1970). In a subsequent clustering step, a binary look-up table is computed by evaluating the misorientation of each pair of sampled orientations. If this misorientation is below a chosen threshold value, the corresponding entry in the look-up table is set to 1 otherwise to 0. To assign orientations to grains, each grain is associated with a number of so-called elementary volumes according to their size, which is used to find orientations in the look-up table with at least this number of orientations having low misorientation. The corresponding crystallographic grain orientation is then the average of orientations with low misorientation to each other. The parameters, i.e., the number of elementary volumes per grain and the threshold value for the misorientation, have to be chosen judiciously. The TOP method is implemented in Python with Cython extension following the optimization procedure outlined in section 4.2.2. Unless otherwise specified, we use a tolerance of tol = 10<sup>−</sup><sup>8</sup> to solve the optimization problem and consider texture coefficients up to rank six. The material model described in Sec. 2.2.2 is implemented in a usermaterial-subroutine (UMAT). The coefficients of the elastic stiffness

tensor are taken from the literature (Wu et al., 2017; Schäfer et al., 2019a), whereas the critical resolved shear stress, assumed to be identical for all slip systems, and the parameters of the kinematic hardening model were fitted to experimental stress-strain hysteresis of the steel C45 using Bayesian optimization (Kuhn et al., 2021). The complete set of used model parameters is summarized in Tab. 4.1.


**Table 4.1:** Parameters used for the crystal plasticity model (Wu et al., 2017; Schäfer et al., 2019a).

To efficiently compute the effective stiffness as well as the macroscopic stress-strain hysteresis, we use the FFT-based solver FeelMath (Fraunhofer ITWM, 2021; Wicht et al., 2020a;b). For the stiffness computations we rely on the conjugate gradient method (Zeman et al., 2010; Brisard and Dormieux, 2010), whereas for the non-linear problem we use a Newton-CG method (Gélébart and Mondon-Cancel, 2013; Kabel et al., 2014). For both problems we use the Moulinec-Suquet discretization (Moulinec and Suquet, 1994; 1998). For a perspective of solution schemes and discretizations, we refer to the recent review article by Schneider (2021). By default, we carry out the computations on periodic microstructures, discretized by 64<sup>3</sup> voxels. Please note that we apply periodic boundary conditions to compute the stiffness and hysteresis.

#### **4.3.2 Linear elastic stiffness**

In this section, we study the effect of different orientation-sampling techniques on the effective stiffness of polycrystalline microstructures. In order to minimize the influence of the underlying microstructure morphology, we use a fixed grain microstructure for each realization and *all* orientation sampling methods. This is illustrated in Fig. 4.2, where we show the results of different sampling techniques for a fixed grain structure with grains of identical volume.

**Figure 4.2:** Furnishing a grain microstructure with 128 grains of identical volume with orientations sampled by three different methods. The colors correspond to the inversepole-figure color key in 010 direction (Nolze and Hielscher, 2016). (a) TOP, (b) Neper, (c) Random.

**Uniform CODF** We start with the case of a unique grain-size and a uniform orientation distribution, corresponding to mechanically isotropic behavior (Krawietz, 1999; Bertram et al., 2000; Böhlke and Bertram, 2001). For the results to be representative, it is necessary to determine the number of grains which ensure an isotropic effective material response, see for example Kanit et al. (2003) and Yang et al. (2019). In this spirit, we investigate microstructures with an increasing number of grains and study their effective stiffness.

As discussed in Sec. 4.2.1, for a uniform CODF, all texture coefficients vanish, i.e.,

$$
\bar{\mathbb{V}}'\_{\langle \beta \rangle} = 0 \tag{4.18}
$$

holds for all considered texture coefficients. To quantify the anisotropy of the stiffness tensor we compare to the best approximation by an isotropic tensor (see equation (4.22)), i.e., we project the computed stiffness tensor onto the space of isotropic tensors of fourth order. For a detailed discussion see the work by Fedorov (1968) and Arts (1993). We compute the mean stiffness

$$\overline{\mathbb{C}}\_{G} = \frac{1}{N} \sum\_{n=1}^{N} \mathbb{C}\_{G,n} \tag{4.19}$$

of = 10 realizations and extract the Lamé constants via

$$
\mu^{\rm app} = \frac{1}{3} \left( \bar{C}\_{G,44} + \bar{C}\_{G,55} + \bar{C}\_{G,66} \right) \tag{4.20}
$$

$$\lambda^{\rm app} = \frac{1}{6} \left( \bar{C}\_{G,12} + \bar{C}\_{G,13} + \bar{C}\_{G,23} + \bar{C}\_{G,21} + \bar{C}\_{G,31} + \bar{C}\_{G,32} \right), \tag{4.21}$$

where ¯*,* denotes the -th entry of the stiffness tensor in Voigt notation (Cavallini, 1999).

Using the best isotropic approximation C iso ( app*,* app) based on the extracted Lamé constants, we introduce the isotropy error iso via

$$\delta^{\rm iso} \left( \mathbb{C}\_{G,1}, \dots, \mathbb{C}\_{G,N} \right) = \frac{\left\| \mathbb{C}^{\rm iso} \left( \mu^{\rm app}, \lambda^{\rm app} \right) - \overline{\mathbb{C}}\_{G} \right\|}{\left\| \overline{\mathbb{C}}\_{G} \right\|}, \tag{4.22}$$

measuring the degree of anisotropy present in the computed stiffness. For an increasing number of grains

$$G \in \{32, 64, 96, 128, 256, 384, 512, 768, 1024, 1536\},$$

we show the resulting isotropy error for the three different orientation sampling methods in Fig. 4.3a. We observe a decreasing isotropy error for all methods with an increasing number of considered grains. All methods decrease the isotropy error at a similar rate. However, they differ in the initial error level. For instance, all sampling techniques reach a low isotropy error for 1536 grains, namely 0*.*251%, 0*.*041% and 0*.*027% for random sampling, the Neper and TOP methods, respectively. To reach a mean error below 1%, the microstructure has to consist of more than 64 grains if the orientations are sampled randomly or generated by Neper. For all investigated grain counts, the TOP method produces the lowest isotropy error. Neper starts with a substantially higher error (by roughly one order of magnitude) at low grain counts and reaches a similar performance to TOP for more than 300 grains. For the naive random sampling, the isotropy error has a quite large offset to the more involved algorithms.

**Figure 4.3:** Isotropy and total error for effective stiffnesses computed from microstructures with uniform orientations and unique grain size distribution. (a) Isotropy error, (b) Total error.

In addition to evaluating the degree of isotropy, we investigate the deviation from the effective, infinite-volume stiffness. As our ground truth, we consider the mean of ten apparent stiffnesses, each computed using volume elements consisting of 10 000 grains and discretized by 128<sup>3</sup> voxels (see Fig. 4.4a).

**Figure 4.4:** One realization of microstructures with with TOP based orientations. The color corresponds to the ipf color key in 100 direction. (a) 10 000 grains, (b) 1024 grains.

The orientations are sampled using the TOP method. The mean stiffness and the 95% confidence intervals, computed via Student's t-distribution (Student, 1908), are given in Tab. 4.2. The isotropy error of the mean stiffness is iso = 0*.*012%.


**Table 4.2:** Mean and 95% confidence intervals in GPa for the stiffness in Voigt'snotation computed by averaging ten realizations of microstructures with 10 000 grains and uniformly distributed TOP orientations.

We define the total error tot as the mean relative error between the stiffness of each realization and the one given in Tab. 4.2, i.e.,

$$\delta^{\text{tot}} = \frac{1}{N} \sum\_{n=1}^{N} \frac{||\mathbb{C} - \mathbb{C}\_{G,n}||}{||\mathbb{C}||}. \tag{4.23}$$

For the total error, shown in Fig. 4.3b, we make similar observations as for the isotropy error. All of the methods decrease the total error at a similar rate, but differ initially. For the case of orientations generated by Neper, the error for 32 grains is tot = 1*.*17% and therefore roughly twice as large compared to the TOP method with tot = 0*.*521%. To reach a similar error with randomly sampled orientations about 768 grains have to be considered. Random sampling leads to a mean error of 0*.*677% for 1536 grains. The errors produced by the Neper and TOP method are similar to each other with 0*.*073% and 0*.*078%, respectively.

To understand the similar rates of error decrease more thoroughly, it is helpful to decompose the total error tot into two contributions (Kanit et al., 2003). The first part is the random error and quantifies the inaccuracy associated with working on a reduced representation of the ground truth. The second contribution quantifies artificial long-range correlations introduced by working on periodic microstructures (Kanit et al., 2003; Schneider et al., 2021). We attribute the visible offset in Figs. 4.3 and 4.6 to the random error, as we use the same geometric representations for each orientation sampling method. Thus, a smaller random error is achieved by the TOP method and further reduction of the total error tot is attributed to increasing the cell-size, i.e., increasing the number of grains.

Up to this point, our investigations were based on a polycrystal with a unique GSD, i.e., a unique grain size. However, to account for the influence of the grain size on the mechanical response, it may often be necessary, and therefore desirable, to match more realistic grain size distributions when generating synthetic polycrystalline microstructures. Thus, we turn to polycrystals with a log-normal GSD, as typically observed in real-world samples (Spettl et al., 2014), with a mean equivalent diameter equal to unity and a standard deviation of 0*.*15. Fig. 4.5 shows an example of a microstructure consisting of 128 grains, equipped with orientations from the three different sampling techniques.

**Figure 4.5:** Polycrystalline microstructure realizations with 128 grains following a lognormal grain size distribution and orientations generated by three different methods. The colors correspond to the ipf color key in 010 direction. (a) TOP, (b) Neper, (c) Random.

The isotropy as well as the total error for different sampling methods are shown in Fig. 4.6. For the log-normal GSD, we register a notable decrease in accuracy for the Neper sampling method compared to the unique grain size. Considering the case of 32 grains, the isotropy error and total error increase from 1*.*22% and 1*.*17% to 1*.*62% and 2*.*19%, respectively. This loss of accuracy persists, even for larger grain counts, e.g., for 1536 the total error for the unique and log-normal GSD are 0*.*07% and 0*.*37%, respectively. For randomly sampling the Haar distribution, the influence of a log-normal grain size distribution is smaller. For instance the biggest difference in the isotropy error is 1*.*23% and 1*.*60% for the unique and log-normal GSD and 32 grains, respectively. The proposed TOP method takes the volume fraction of each grain into account in an explicit way when optimizing the orientations. This results in strikingly similar error levels for both the unique and the log-normal case. Whereas the total error values realized by microstructures with 32 grains differ slightly for the unique and log-normal case, the resulting isotropy error is iso = 0*.*16% for both GSDs.

**Figure 4.6:** Isotropy and total error for the effective stiffnesses computed from microstructures with uniform orientations and log-normal grain size distribution. (a) Isotropy error, (b) Total error.

We investigate the influence of the maximum rank of the texture coefficients considered in our optimization scheme in the case of a uniform orientation distribution. For this purpose, we consider the case of ten microstructures consisting of 1024 grains, see Fig. 4.4b for an example of one realization. We study two different cases: Using solely the texture coefficient of rank four to optimize the orientations and considering the coefficients of rank four *and* six. For these two cases, the isotropy and total error are shown in Fig. 4.7 for different values of the tolerance tol used in the optimization procedure. Both errors show a decrease up to a value of tol = 10<sup>−</sup><sup>6</sup> , after which the resulting errors do not change. Interestingly, considering only the texture coefficient of rank 4 appears to be beneficial. As all total errors are below 0*.*1% and all isotropy errors even below 0*.*03%, using texture coefficients of rank four and solving the problem up to a tolerance of tol = 10<sup>−</sup><sup>6</sup> is sufficient for the case of a uniform orientation distribution when considering 1024 grains and solely the macroscopic stiffness is of interest.

**Figure 4.7:** Comparison of the tot and iso for the stiffness computed for 1024 grains with TOP orientations and different considered texture coefficients plotted on a finely resolved -axis. (a) Isotropy error, (b) Total error.

**Textured CODF** To further investigate the capabilities of the TOP method, we turn to a non-uniform CODF, i.e., a *textured* polycrystal. The prescribed CODF was generated by MTex (Hielscher and Schaeben, 2008), taken from the MTex documentation (MTex, 2017), see Fig. 4.8 for the corresponding pole figures. As MTex allows the sampling of CODFs, we draw 50 000 samples at random for computing the texture coefficients, assuming the same weight for each sample.

As a ground truth we define the mean stiffness of ten realizations, each with 10 000 grains. The resulting stiffness for TOP orientations is given, with its respective 95% confidence intervals, in Tab. 4.3. The isotropy error of this stiffness computes to iso = 5*.*54%, i.e., a slight anisotropy appears.

**Figure 4.8:** Pole figures of the (generated) textured CODF (MTex, 2017).


**Table 4.3:** Mean and 95% confidence intervals in GPa for the stiffness in Voigt notation computed by averaging ten realizations of microstructures with 10 000 grains and TOP orientations for a synthetic CODF.

For this texture, we investigate the approximation quality of the stiffness for a varying number of grains, each with identical volume. The total error for randomly sampling from the generated orientations and using texture coefficients is shown in Fig. 4.9a. Randomly sampling the generated orientations gives a total error of tot = 4*.*47% for 32 grains and reaches an error below 1% for 936 orientations. The mean total error achieved by the TOP method of tot = 0*.*51% for 32 grains actually lies below the error value achieved by randomly sampling 1536 orientations. For the latter number of grains, TOP achieves an error tot = 0*.*07%. This difference is attributed to the notable offset between the random sampling and TOP method, as both decrease tot with the same rate. Let us consider the case of a log-normal grain size distribution. For the TDT algorithm we assign eight elementary volumes to the smallest grain

and increase the number of elementary volumes for each grain according to its size. We set the threshold misorientation value to 5 ∘ . In Fig. 4.9, we provide the total error tot for ten realizations with a varying number of grains. For the case of randomly sampling from given orientations, we observe a slight increase in the error value induced by the underlying log-normal grain size distribution. For instance, for a microstructures with 32 grains, the mean error is tot = 3*.*12% and tot = 3*.*74% for the unique and log-normal GSD, respectively. This effect decreases when a larger number of grains is considered, as the effect of a single, large grain with specific orientation on the overall response decreases. In contrast, the TOP method is not adversely affected. Indeed, during optimization, the volume fraction is explicitly taken into account when computing the texture coefficients, see equation (4.12). Using the TDT algorithm results in a lower total error than random sampling for all grain numbers considered. For instance, for a 64-grain microstructure, the total error is tot = 3*.*75% and tot = 2*.*77% for random sampling and the TDT, respectively. For both algorithms the error decreases with a similar rate as for the proposed TOP algorithm, whereas they both result in higher total errors than using TOP.

**Figure 4.9:** Total error tot of the effective stiffness computed for microstructures with a unique and log-normal grain size distribution and a synthetic CODF. (a) Unique grain size distribution, (b) Log-normal grain size distribution.

**Highly textured CODF** In practical applications, e.g., cold rolled steel, the intensities in the pole figure may reach values as high as ten. To investigate this scenario, we next consider a case with an increased texture in the CODF. We rely on synthetically generating a CODF using MTex (Hielscher and Schaeben, 2008) and show the resulting pole figures in Fig. 4.10.

**Figure 4.10:** Pole figures of the (generated) CODF with increased texture (MTex, 2017).

For the ground truth we proceed in the same way as for the slightly textured CODF, using ten microstructures with 10 000 grains of equal volume, equipped with orientations from the TOP method to compute the mean stiffness. For this case, the mean and the 95% confidence intervals are given in Tab. 4.4. The isotropy error is iso = 18*.*78% which is more than three times the error of the slightly textured case, i.e., iso = 5*.*54%.


**Table 4.4:** Mean and 95% confidence intervals in GPa for the stiffness in Voigt notation computed by averaging ten realizations of microstructures with 10 000 grains and TOP orientations for a synthetic CODF.

First, we investigate the case of a unique grain size distribution and show the resulting total error in Fig. 4.11a. For TOP and random sampling, the error decreases with a similar rate, which is consistent with our observations in the slightly textured case. For TOP as well as for random sampling the total error is slightly lower than for the previously investigated CODF, e.g., for 32 grains the total error is tot = 3*.*62% and  tot = 0*.*40% for random sampling and TOP, respectively. This holds for higher grain numbers as well. Indeed, for 1536-grain microstructures, randomly sampling orientation data leads to a total error of tot = 0*.*53%, whereas using orientations generated by TOP results in an error of tot = 0*.*06%. However, the relative difference between the total errors for the two CODFs is *lower* for random sampling than for TOP.

**Figure 4.11:** Total error tot of the effective stiffness computed for microstructures with a unique and log-normal grain size distribution and a highly textured CODF. (a) Unique grain size distribution, (b) Log-normal grain size distribution.

For the case of a log-normal GSD, compared to microstructures with grains of equal volume, the total error increases for *both* methods. For microstructures with 64 grains equipped with orientations randomly sampled from experimental data, the total error is tot = 2*.*82% and tot = 3*.*11% for a uniform and log-normal GSD, respectively. This observation holds for the TOP method as well, e.g., using 64 grains leads to an error increase from tot = 0*.*40% for a unique GSD to tot = 0*.*45% if the grain sizes follow a log-normal distribution. For the TDT algorithm and a grain count below 768, we set the number of elementary volumes

for the smallest grain to eight. To account for the increased grain count, we increase the number of elementary volumes to twelve for 1024 and 1536 grains in the microstructure, whereas we retain the threshold of 5 ∘ for the misorientation computations. For our choice of parameters and grain numbers up to 256, we observe that the resulting error is close to random sampling. For instance, the total error obtained using a 256-grain microstructure is tot = 0*.*94% and tot = 1*.*37% for orientations from random sampling and the TDT algorithm, respectively. When increasing the grain count, there seems to be a limiting accuracy that the TDT algorithm can reach. Indeed, the total error does not decrease below 1%. The source of this phenomenon needs to be investigated more thoroughly, and is beyond the scope of this work.

#### **4.3.3 Cyclic stress-strain hysteresis**

We expand our investigation into the elasto-plastic regime, focusing on the effect of the orientation sampling method on the cyclic stress-strain hysteresis of the material. As boundary condition, we use a macroscopic strain which follows a triangular path with an amplitude of = 0*.*7% and a cycle time of four seconds. To ensure a stabilized cyclic stress-strain hysteresis, we compute two cycles in total and use the last one as our quantity of interest (Schäfer et al., 2019a).

Because of the increased computational cost, we restrict the investigations to grain counts

$$G \in \{32, 64, 128, 256, 512, 1024\}$$

and use five realizations per number of grains, i.e., = 5. We use the material parameters specified in Tab. 4.1.

**Uniform CODF** For the uniformly distributed orientations, we asses the isotropy of the results. For this purpose, and a single realization, we compute three load cases to obtain the cyclic stress-strain hysteresis

in three different directions, i.e., -,  - and -direction. For a perfectly isotropic response, the stress values would coincide for every considered direction. Thus, to measure the deviation from this isotropic result, we use the average stress values

$$\bar{\sigma}\_s = \frac{1}{3} (\sigma\_{XX,s} + \sigma\_{YY,s} + \sigma\_{ZZ,s}) \tag{4.24}$$

of all directions at each time step as our reference. Then, for each realization , we compute the sum of the squared relative differences between the stresses in every direction and the mean of all directions weighted by the number of stress values , i.e.,

$$\delta\_n^{\text{hys,iso}} = \sqrt{\frac{1}{S} \sum\_{s=1}^{S} \left( \frac{\sigma\_{XX,s}}{\bar{\sigma}\_s} - 1 \right)^2 + \left( \frac{\sigma\_{YY,s}}{\bar{\sigma}\_s} - 1 \right)^2 + \left( \frac{\sigma\_{ZZ,s}}{\bar{\sigma}\_s} - 1 \right)^2}. \tag{4.25}$$

Equation (4.25) measures the mean relative deviation in all directions from the mean stress value (4.24). The quantity hys,iso extends the isotropy error defined for the stiffness, see equation (4.22), where the ideal isotropic case corresponds to the mean stress values in all directions. We compute the mean error of all realizations = 5 by

$$\delta^{\text{hys,iso}} = \frac{1}{N} \sum\_{n=1}^{N} \delta\_n^{\text{hys,iso}} \tag{4.26}$$

to get confidence in our results.

We use the microstructures from Sec. 4.3.2 with orientations prescribed by Neper, TOP and random sampling. For an increasing number of grains, we show the error hys,iso in Fig. 4.12a. Neper and TOP behave similar to each other, both lying below the error values achieved by random sampling. For example, the mean errors obtained from microstructures with 32 grains are hys,iso = 5*.*63%, hys,iso = 5*.*97% and

 hys,iso = 1*.*10% for random sampling, Neper and TOP, respectively. For an increasing number of grains in the microstructure and randomly sampled orientations, the error decreases more slowly than for the other methods. Also, random sampling results in the highest mean error value of hys,iso = 4*.*49% for 1024 grains. We observe a steeper decrease in hys,iso for an increasing number of grains for Neper and TOP, both lying close to each other. For example, microstructures with 64 grains produce an error of hys,iso = 1*.*19% and hys,iso = 1*.*40% for TOP and Neper orientations, respectively. For 1024 grains, the error levels are hys,iso = 0*.*27% and hys,iso = 0*.*16% for Neper and TOP orientations, respectively.

**Figure 4.12:** Hysteresis isotropy error, equation (4.25) for a varying number of grains . (a) Different orientation generation methods, (b) Different considered texture coefficients for TOP.

Following the procedure in Sec. 4.3.2, in addition to investigating the degree of isotropy, we would like to assess the ability to reproduce the effective mechanical response with a minimum number of grains. In the case of non-linear mechanical behavior, we define our ground truth as the stress-strain hysteresis computed for five realizations of

a microstructure with 10 000 grains, discretized by 128<sup>3</sup> voxels. From these five realizations, we compute the mean stress-strain curve as

$$
\hat{\sigma}\_{r,s} = \frac{1}{N} \sum\_{n=1}^{N} \sigma\_{r,s,n},
\tag{4.27}
$$

where *,,* denotes the macroscopic stress value in direction at a given time step and refers to the number of realizations, i.e., = 5. For the considered loading directions, the resulting stress-strain hysteresis are shown in Fig. 4.13. We observe that the individual curves lie on top of each other, i.e., there is no anisotropy present.

**Figure 4.13:** Comparison of mean macroscopic stress-strain hysteresis in different loading directions computed using five microstructures with 10 000 grains and uniformly distributed TOP orientations. (a) and  , (b) and , (c)  and .

For each number of grains and each realization, we compute the root of the mean squared relative error in each direction as

$$\delta\_{n,r}^{\text{hys,gt}} = \sqrt{\frac{1}{S} \sum\_{s=1}^{S} \left( \frac{\sigma\_{s,r,n}}{\hat{\sigma}\_{s,r}} - 1 \right)^2},\tag{4.28}$$

where refers to the loading direction, i.e., ∈ {*, ,* } in our case. Then, we compute the mean error over all considered directions

$$\delta^{\text{hys,gt}} = \frac{1}{N} \sum\_{n=1}^{N} \frac{1}{R} \sum\_{r=1}^{R} \delta\_{n,r}^{\text{hys,gt}},\tag{4.29}$$

where denotes the number of considered directions, i.e., = 3 in this case.

Comparing hys,gt for different orientation sampling methods in Fig. 4.14, we observe similar trends as for the isotropy error hys,iso. For all methods, the error decreases with an increasing number of grains in the microstructure. For a small number of grains, the error resulting from TOP orientations is smallest with hys,gt = 1*.*80% and hys,gt = 0*.*78% for 32 and 256 grains, respectively. The error from using Neper orientations is higher, with hys,gt = 3*.*93% and hys,gt = 0*.*82%. Randomly sampling the Haar distribution results in an error of 9*.*00% and 2*.*15% for the same number of grains. The error for all three methods and 1024 grains are 0*.*43%, 0*.*44% and 1*.*18% for TOP, Neper orientations and random sampling, respectively.

The error hys,gt for taking only the texture coefficient of rank four into account, is shown in Fig. 4.14b together with the previously discussed results for considering texture coefficients with rank four and six. We observe that the error when accounting solely for rank four texture coefficients is higher than the error produced when considering higher ranks. The difference is less pronounced than for hys,iso .

**Figure 4.14:** Hysteresis total error, equation (4.29) for a varying number of grains with uniform CODF and unique GSD. (a) Different orientation generation methods, (b) Different considered texture coefficients for TOP.

To extend our studies to a non-unique grain size distribution, we use microstructures with a log-normal grain size distribution. We fix the mean and standard deviation to mean = 1 and stdev = 0*.*15, respectively. To reduce the computational effort and because Neper and TOP provided the most promising results, we only consider orientations generated by Neper and TOP in the following.

Fig. 4.15a shows hys,iso for an increasing number of grains . For every number of grains, the TOP methods provides a smaller error compared to Neper. For instance, using 32 grains, the error for TOP and Neper is hys,iso = 3*.*72% and hys,iso = 12*.*02%, respectively. The influence of the underlying GSD manifests. Indeed, for both cases, the values are larger than for the unique grain size distribution.

Similar to the uniform GSD, the hysteresis error closely follows the trend observed fo the isotropy error, see Fig. 4.15b.

**Figure 4.15:** Comparison of hys,gt and hys,iso for varying number of grains with uniform CODF and log-normal GSD. (a) Hysteresis isotropy error, equation (4.25), (b) Hysteresis total error, equation (4.29)

**Textured CODF** For the case with mild anisotropy, we consider the synthetic CODF described in Sec. 4.3.2. We compute the stress-strain hysteresis using five microstructures consisting of 10 000 grains discretized by 128<sup>3</sup> voxels, see Fig. 4.4a. In accordance with the case of a uniform orientation distribution, we use the mean stress values of these five realizations as our ground truth. The resulting stress-strain hysteresis are shown in Fig. 4.16 for all three considered loading directions. We observe a slight anisotropy in  -direction, whereas the stress-strain curves in - and -direction coincide.

**Figure 4.16:** Comparison of mean macroscopic stress-strain hysteresis in different directions computed for a microstructure with 10 000 grains and a slightly textured CODF. (a) and  , (b) and , (c)  and .

We show the total error to the mean stress values, i.e., equation (4.29), in Fig. 4.17a for the TOP method as well as for randomly sampling from given orientations.

Using random orientation sampling produces a larger error for all grain numbers considered. Especially for a small number of grains, the TOP method results in a visibly smaller error than for random sampling the experimental data. For 32 grains, the hysteresis total error hys, tot is 3*.*60% and 1*.*28% for random sampling and TOP, respectively. For random sampling, the error reduces to 0*.*76% for 1024 grains, which is close to the value achieved by TOP, with hys,gt = 0*.*52%.

**Figure 4.17:** hys,gt for varying number of grains with a slightly textured CODF. (a) Unique GSD, (b) Log-normal GSD.

We observe similar behavior for the case of a log-normal grain size distribution with mean = 1 and stdev = 0*.*15 in Fig. 4.17b. The error for 32 grains increases for both kinds of orientations sampling methods, namely to 3*.*84% and 1*.*41% for random and TOP sampling, respectively. For both the unique and log-normal case, similar errors of 0*.*85% and 0*.*44% are achieved for randomly sampling experimental orientations and using 1024 grains. Interestingly, the error for a log-normal GSD is actually smaller than for the unique GSD. For instance using the TOP method and 512 grains, we observe an error of hys,gt = 0*.*66% and hys,gt = 0*.*49% for the unique and log-normal distributions, respectively. For the TDT algorithm, we observe a lower error than for random sampling when an intermediate number of grains is considered, i.e., for grain counts of 64, 128 and 256. For instance, the error for a microstructure consisting of 64 grains equipped with orientations of the TDT algorithm is hys,gt = 2*.*23% and hys,gt = 2*.*22% for 128 grains. With 64 and 128 grains, the random sampling leads to an error of hys,gt = 3*.*84% and

 hys,gt = 2*.*80%, respectively. For a larger number of grains, above 256, the error computed for the TDT algorithm exceeds the error of the randomly sampled orientations from given data. For 1024 grains, the error for random sampling is hys,gt = 0*.*85% and hys,gt = 1*.*00% for the TDT algorithm. All of the observed error values are above the errors obtained by TOP, e.g., using 32 grains, the error for the TOP method is hys,gt = 1*.*42% whereas the TDT algorithm and random sampling lead to an error of hys,gt = 4*.*37% and hys,gt = 3*.*84%, respectively.

**Highly textured CODF** Last but not least, we consider the GSD outlined in Sec. 4.3.2 with an increased degree of anisotropy. Similar to the case of a slight anisotropy, we compute the stress-strain hysteresis for five microstructures consisting of 10 000 grains and a discretization of 128<sup>3</sup> voxels. As our ground truth we use the mean stress of these five realizations.

For microstructures with grains having a unique grain size distribution and orientations from TOP or randomly sampling experimental data, we show the total error in Fig. 4.18a. We observe that, for all grain counts considered, using the TOP method results in lower error values compared to randomly sampling from given orientation data. For instance, using microstructures with 32 grains leads to a total error of hys,gt = 7*.*54% and hys,gt = 17*.*44% for TOP and random sampling, respectively. Thus, we observe an increase in the total error in comparison to the slightly textured CODF for both sampling methods and all microstructures. Indeed, for a 1024-grain microstructure equipped with orientations from TOP, the error increases from hys,gt = 1*.*42% for the slightly textured case to hys,gt = 5*.*24% for the case of higher texture. This observation holds for randomly selecting orientations from given orientation data, e.g., for a microstructure consisting of 256 grains the error for the highly textured CODF is hys,gt = 10*.*1% whereas it is hys,gt = 1*.*82% for the slightly textured case.

**Figure 4.18:** hys,gt for varying number of grains with a highly textured CODF. (a) Unique GSD, (b) Log-normal GSD.

Fig. 4.18b shows the total error for the case of a log-normal GSD equipped with orientations from TOP, TDT and random sampling. We make similar observations to the slightly textured CODF, i.e., an increase in the total error compared to the results for the unique GSD. For instance, a 64-grain microstructure with TOP orientations leads to an increase in the total error from hys,gt = 5*.*36% for a unique GSD to hys,gt = 9*.*62% for a log-normal GSD. Comparing the same microstructures equipped with orientations from randomly sampled orientation data, the error increases from hys,gt = 13*.*2% to hys,gt = 16*.*3% for a unique and a log-normal GSD, respectively. For lower grain counts, i.e., below 256, using the TDT algorithm results in similar error values as randomly sampling orientation data. An exception is the 32-grain microstructure, for which the total error value hys,gt = 11*.*7% is close to the value obtained using TOP, i.e., hys,gt = 8*.*02%. For higher grain counts, i.e., above 256, the TDT algorithm does not decrease the total error but instead we observe an increase in the error values obtained, in line with

observations made investigating the performance for the linear elastic properties in Sec. 4.3.2.

## **4.4 Conclusion**

In this work, we proposed to use the coefficients of a tensorial Fourier expansion of the crystallite orientation distribution function (Guidi et al., 1970) to equip digital polycrystalline microstructures with crystallographic orientations for micromechanical simulations. Our proposed method is based on minimizing the difference between the current and the prescribed tensorial Fourier, or texture, coefficients and uses a gradient descent scheme on the Lie group (3).

We compared the proposed texture optimization for prescribing orientations (TOP) method to different state-of-the-art methods, e.g., implemented in the sophisticated microstructure generation tool Neper (Quey et al., 2011; 2018). In a first step, we investigated the homogenized stiffnesses of polycrystals for the case of a uniform and two textured crystallite orientation distribution functions (CODF). Subsequently, we extended our studies to the non-linear case, where we investigated the macroscopic cyclic stress-strain hysteresis of the microstructures. For both, the linear elastic and the non-linear plastic case, we considered a unique as well as a log-normal grain size distribution (GSD). By introducing suitable error measures we investigated and compared the performance of the proposed method.

In the isotropic, linear elastic case, TOP provided better results compared to the Neper method and random sampling of orientations. Using TOP, an isotropic effective stiffness could already be achieved for small grain numbers. Owing to the fact that the volume fraction of each individual grain is explicitly accounted for, the advantage becomes more pronounced when dealing with microstructures having a log-normal grain size distribution. Additionally, with TOP it is possible to reproduce the

linear-elastic behavior of polycrystals with a unique GSD and a textured CODF more accurately and efficiently than via a random sampling from experimental orientation data. This holds as well for the case of a lognormal grain size distribution. Comparing to the Texture Discretization Technique (TDT) algorithm proposed by Melchior and Delannay (2006), which also considers grain size during orientation assignment, the TOP method performed better for both CODFs considered. Our intensive numerical studies revealed that the performance of the TDT algorithm, in our setting, critically depends on the choice of parameters, i.e., the misorientation value and the number of elementary volumes per grain. For the non-linear plastic behavior, the results of the Neper method were very similar to the ones provided by TOP, showing the capabilities of the dedicated algorithm. Although the effect was less pronounced than for the case of linear-elastic behavior, we observed that a underlying lognormal GSD results in a decreased performance for the Neper method. TOP, on the other hand, was able to produce similar results as for the unique GSD. In addition, microstructures with orientations provided by the TOP method allow to accurately compute the effective, nonlinear behavior of polycrystals with an underlying texture. For both textured CODFs considered in this work, the defined error measures were lower compared to randomly sampling orientation data. The errors for a log-normal GSD obtained by TOP were below the ones obtained with random sampling or the algorithm of Melchior and Delannay (2006), even for a small number of grains. Investigating a highly textured CODF and physically non-linear visco-plastic behavior, we observed higher error values compared to the slightly textured CODF for all algorithms considered. This contrasts with the linear elastic case where we observed a smaller error for the highly textured CODF. However, using orientations generated by the TOP method leads to significantly smaller errors than using the other two algorithms. Concerning the maximum number of texture coefficients which should be taken into consideration, we observed that, for the linear elastic case, a low tolerance and only the texture coefficient of rank four a sufficient. For the non-linear behavior, we observed that accounting for the texture coefficient of rank six is beneficial.

As the computational effort of micromechanical studies is mainly dominated by computing the effective behavior we omitted a comparison of the computational performance of TOP to the other methods.

To conclude, we showed that extracting relevant data from the CODF in terms of tensorial texture coefficients leads to the most flexible and performing method for generating crystallite orientations for digital representations of polycrystalline microstructures. As orientation assignment is typically treated as a post-processing step in microstructure generation, it is possible to couple the proposed algorithm with well-established microstructure generators (Groeber and Jackson, 2014; Prasad et al., 2019; Henrich et al., 2020). With this modular structure, it is possible to use TOP for generating polycrystalline representations for a variety of applications (Flipon et al., 2020b; Vajragupta et al., 2020).

As an additional benefit, for generating orientations, the TOP method requires only nine (or 22 variables, depending on the highest texture coefficient rank considered) to be stored, contrasting with methods that rely on the entire experimental database. Because of this low number of parameters, it is possible to fuel data driven methods (Pütz et al., 2020; Gajek et al., 2021). Additionally, as experimental data is always afflicted with some degree of measurement uncertainty, investigating the influence of the texture on the overall macroscopic response might be an interesting topic, i.e., via uncertainty quantification (Bandyopadhyay et al., 2019; Kasemer et al., 2020).

### **Chapter 5**

# **Identifying material parameters in crystal plasticity by Bayesian optimization <sup>1</sup>**

## **5.1 Introduction**

Using suitable methods to generate representations of polycrystalline microstructures, see Chap. 3 and Chap. 4, we still need a compliment material model to predict the cyclic behavior of polycrystalline metals. The crystal plasticity method, outlined in Sec. 2.2.2, proves to be a powerful method but needs suitable model parameters to capture the relevant features of the cyclic deformation behavior.

Within this chapter we propose using Bayesian optimization with Gaussian processes for calibrating single crystal parameters inversely based on polycrystalline experiments. We improve upon related research (Schäfer et al., 2019a) by assessing the capabilities of our proposed framework by comparing it to a number of extremely powerful methods, including a genetic algorithm and modern derivative-free optimization schemes. This chapter is structured as follows. We recall Gaussian-process based Bayesian optimization in section 5.2.

<sup>1</sup> This chapter is based on the publication by Kuhn et al. (2021) whereas minor typographical and formatting changes have been made for cohesion of the manuscript and the readers convenience

We investigate the performance of the proposed approach in section 5.3, carefully studying the optimal algorithmic parameters, comparing to state-of-the-art genetic algorithms (Schäfer et al., 2019a) as well as powerful derivative-free optimization methods (Rios and Sahinidis, 2013; Huyer and Neumaier, 1999; 2008) and demonstrating the applicability of the method to large-scale polycrystalline microstructures in section 5.3.6.

## **5.2 Bayesian optimization**

Bayesian optimization (BO) is an approach for solving optimization problems where the objective function is expensive to compute and the gradient of the objective function is not available. Such optimization problems occur, for instance, when evaluating the function corresponds to running a (possibly large-scale) simulation, as is common in virtual crash tests (Raponi et al., 2019) or for complex chemical reactions (Häse et al., 2018). Such problems might also be approached by automatic differentiation (Griewank and Walther, 2008), which might prove difficult to integrate into an existing simulation code, or by evolutionary algorithms, like particle swarm algorithms (Blum and Merkle, 2008). However, the latter are typically limited to small spatial dimensions and inexpensive function evaluations.

In contrast, Bayesian optimization constructs a surrogate model of the optimization problem based on probabilistic ideas, accounting for uncertainty by Bayesian statistics. We restrict to GAUSSIAN PROCESS REGRESSION as our uncertainty model. For other approaches we refer to Shah et al. (2014) and Kushner (1964).

Suppose that we are concerned with the optimization problem

$$f(x) \longrightarrow \min\_{x \in A},\tag{5.1}$$

where ⊆ R is a non-empty set in spatial dimensions whose membership is easily tested (for instance, a box). For (Gaussian process based) Bayesian optimization, we need to specify:


$$
\mu: A \to \mathbb{R} \quad \text{and} \quad \sigma: A \to \mathbb{R}\_{\geq 0},
$$

which we will refer to as (the local estimates of) the mean and standard deviations of the function values at each ∈ .

There are specific prerequisites for sensible kernel functions. We refer to section 4 in Rasmussen-Williams Williams and Rasmussen (2006) for details, and will discuss our choice below.

For given , and , Bayesian optimization proceeds as follows. Suppose that the objective function was evaluated at points 1*, . . . ,*  in already, i.e., the values (1)*, . . . ,* () are known. In the Bayesian approach, it is assumed that the vector

$$f\_n = \begin{bmatrix} f(x\_1), \dots, f(x\_n) \end{bmatrix}^T \in \mathbb{R}^n$$

was drawn at random from a probability distribution. In Gaussian process regression (Mockus, 1994), the latter vector is assumed to follow a multivariate normal distribution with mean ∈ R and Σ ∈ R ×

given by

$$\mu\_n = \left[ m(x\_1), \dots, m(x\_n) \right]^T \quad \text{and}$$

$$\Sigma\_n = \begin{bmatrix} K(x\_1, x\_1) & K(x\_1, x\_2) & \dots & K(x\_1, x\_n) \\ K(x\_2, x\_1) & K(x\_2, x\_2) & \dots & K(x\_2, x\_n) \\ \vdots & \vdots & \ddots & \vdots \\ K(x\_n, x\_1) & K(x\_n, x\_2) & \dots & K(x\_n, x\_n) \end{bmatrix}. \tag{5.2}$$

To infer the value () at some new point we assume that the joint distribution of the values of over (1*, . . . , ,* ) is also governed by a multivariate normal distribution with mean +1 and Σ+1 for +1 = , as prescribed in equation (5.2). Using Bayes' rule, the conditional distribution of () is also Gaussian with mean

$$
\mu(x) = m(x) + K\_n(x)^T \Sigma\_n^{-1} (f\_n - \mu\_n) \tag{5.3}
$$

and variance

$$
\sigma(x)^2 = K(x,x) + K\_n(x)^T \Sigma\_n^{-1} K\_n(x),
\tag{5.4}
$$

where

$$K\_n: A \to \mathbb{R}^d, \quad x \mapsto \left[ K(x, x\_1), K(x, x\_2), \dots, K(x, x\_n) \right]^T,$$

see section 2 in Williams and Rasmussen (2006). Provided Σ is nondegenerate, and as () corresponds to the -th row of Σ for = 1*, . . . ,* , we immediately see that

$$
\mu(x\_i) = f(x\_i) \quad \text{and} \quad \sigma(x\_i) = 0
$$

hold, i.e., the function interpolates the known values (), and the corresponding standard deviation () vanishes.

Furthermore, () serves as an estimate for the function value () at any potentially newly sampled point ∈ .

In this work, we set the mean function identical to zero, () ≡ 0. As kernel function we consider the anisotropic "Matern5/2" kernel (Matérn, 2013)

$$K\left(x,x\right) = \sigma\_k^2 \left(1 + \sqrt{5}\,\delta\left(x,x\right) + \frac{5}{3}\,\delta\left(x,x\right)^2\right) \exp\left(-\sqrt{5}\,\delta\left(x,x\right)\right) \tag{5.5}$$

in terms of a positive parameter 2 and an anisotropic distance function

$$\delta\left(x\_i, x\_j\right) = \sqrt{\sum\_{d=1}^{D} \frac{(x\_{i,d} - x\_{j,d})^2}{\ell\_d^2}} \quad \text{with} \quad i, j = 1 \ldots n + 1 \tag{5.6}$$

involving dilation parameters *ℓ* for each of the dimensions = 1 *. . .*  of the vector . Indeed, we cannot justify using the squared exponential kernel (Stein, 2012) due to the lack of sufficient smoothness of our objective function . For an overview of alternative kernel functions, we refer to chapter 4 in Williams and Rasmussen (2006). The dilation parameters of the distance function may be interpreted as prefactors used for a nondimensionalization of the variables. To determine the dilation parameters and 2 , a maximum likelihood estimation may be used, see chapter 6 in Bishop (2006). For each new observation the Gaussian process regression is updated, i.e., the parameters of the kernel function (5.5) are determined with respect to *all* previously made observations, continuously improving the model of .

We restricted to discussing Gaussian processes for our noiseless application of Bayesian optimization and refer to Williams and Rasmussen (2006) and Bishop (2006) for a more general discussion.

As the next step, BO searches for an improved guess for the solution of the optimization problem (5.1). A possible approach would be to ignore the estimated statistics and to minimize , exploiting the currently estimated objective values. However, this approach disregards the uncertainty. Indeed, it might happen that the optimum \* is located in regions of high uncertainty. Thus, it might be better to explore first. Acquisition functions provide a suitable exploration-exploitation tradeoff by combining the currently known means and variances into a single function to be optimized.

Denote by \* the lowest objective value observed so far

$$f\_n^\* = \min\_{i=1}^n f(x\_i).$$

Plugging \* and the functions and above into the acquisition function gives rise to the surrogate optimization problem

$$u\_{f\_n^\*,\mu,\sigma}(x) \longrightarrow \max\_{x \in A} \text{.} $$

In contrast to the original optimization problem (5.1), the acquisition function is cheap to evaluate and gradient information is available. Also, the acquisition function provides a trade-off between exploration, i.e., decreasing the uncertainty, and exploitation, i.e., searching in the vicinity of sites with small objective values. Thus, the next point to investigate is selected by maximizing the acquisition function

$$x\_{n+1} = \arg\max\_{x \in A} u\_{f\_n^\*, \mu, \sigma}.$$

Recall that confidence intervals of a normally distributed random variable with mean and standard deviation have the form

$$\left[\mu-\xi\,\sigma,\mu+\xi\,\sigma\right],$$

where is a parameter which determines the probability that a measurement lies in this confidence interval. For instance, a two-sided confidence interval with 95% probability is obtained for ≈ 1*.*96.

For the lower confidence-bound acquisition-function () (Cox and John, 1992), the lower bound of the confidence interval is chosen as the proxy for the objective function , i.e.,

$$u\_{LCB}(x) = -\mu(x) + \xi \,\sigma(x) \tag{5.7}$$

is considered with () and () given by equations (5.3) and (5.4), respectively. This acquisition function tries to improve the currently known least lower bound on .

The quantity may be chosen as a tuneable parameter. Indeed, for small , the acquisition function tends to select points with low mean (). In contrast, high values of encourage the algorithm to choose points where the variance, i.e., the uncertainty, is high. We restrict to a fixed value for and refer to Srinivas et al. (2010) for an adaptive choice of . An alternative strategy, proposed by Mockus et al. (1978) and made popular by Jones et al. (1998), considers the improvement

$$\mathcal{Z}\_n: \mathbb{R} \to \mathbb{R}, \quad f \mapsto \left< f\_n^\* - f \right>\_+,\tag{5.8}$$

where ⟨·⟩<sup>+</sup> = max(0*,* ·) denotes the Macaulay bracket. If ≥ \* , there will be no improvement. For lower values of , the improvement ℐ() increases linearly. Taking the expectation of the improvement (5.8) w.r.t. the normal distribution determined from () and () gives rise to the expected improvement (), which may be expressed analytically in the form

$$u\_{EI}(x) = \left<\Delta\_n(x)\right>\_+ + \sigma(x)\phi\left(\frac{\Delta\_n(x)}{\sigma(x)}\right) - \left|\Delta\_n(x)\right|\Phi\left(\frac{\Delta\_n(x)}{\sigma(x)}\right), \tag{5.9}$$

where Δ() = \* −() and as well as Φ denote the standard normal density and distribution function with mean () and standard deviation (), see Jones et al. (1998).

Lizotte (2008) proposed a shift of equation (5.9)

$$
\Delta\_n(x) = f\_n^\* - \mu(x) - \xi \tag{5.10}
$$

by a parameter . This enables controlling the trade-off between exploitation and exploration and is similar to the parameter used in , see equation (5.7).

The discussed acquisition functions are inexpensive to evaluate, and gradient data is available. Thus, in addition to zeroth order optimization methods, for instance Monte Carlo methods (Wilson et al., 2018) or DIRECT (Brochu et al., 2010), gradient-based solvers like BFGS (Frazier and Wang, 2016; Wang et al., 2018; Picheny et al., 2016) may be used.

Still, finding the global optimum of the acquisition function may be non-trivial. Indeed, even if the original function was convex, the surrogate problem need not be concave. For instance, the expected improvement acquisition function (5.9) has local maxima at all previously explored points 1*, . . . ,* .

Our workflow for maximizing the acquisition function works as follows. First, we sample the acquisition function \* *,,* at 300 points. As our domain of interest is a box, we rely upon the first 300 points of the corresponding Sobol sequence (Sobol, 1967). Then, we select those ten points with the highest acquisition function values, and run BFGS with the selected point as initial point. Finally, we select the maximum value obtained during those BFGS runs.

The employed Bayesian optimization workflow is summarized in Fig. 5.1. Please note that we rely upon a preliminary Latin-hypercube sampling (McKay et al., 2000) to sample the first ten points within the domain of interest .

**Figure 5.1:** Bayesian optimization flowchart.

## **5.3 Computational investigations**

#### **5.3.1 Setup**

We wish to identify material parameters for the quenched and tempered high-strength steel 50CrMo4, as investigated by Schäfer et al. (2019a). The material features a martensitic microstructure and a hardness of 39 HRC. To determine the crystal plasticity parameters of this material, macroscopic strain-driven low-cycle fatigue (LCF) experiments at different strain amplitudes were performed, see Schäfer et al. (2019a) for details on the experimental procedure.

The constitutive material model described in Section 2.2.2 was implemented into a user defined material subroutine (UMAT) within the commercial finite element solver Dassault Systèmes (2018). To thoroughly test the optimization framework, we restrict to the simple microstructure shown in Fig. 5.2, consisting of 8 × 8 × 8 grains, each discretized by 2 3 elements. The different colors in Fig. 5.2 represent distinct, randomly

sampled, orientations. Periodic boundary conditions were used, following the procedure proposed by Smit et al. (1998).

**Figure 5.2:** Simplified volume element used to compute the macroscopic stress strain hysteresis.

The microscopic material parameters we wish to identify are the following:


The remaining model parameters are taken from the literature according to Schäfer et al. (2019a), see Tab. 5.1. Please note that we restrict to slip systems in the lath-plane, as these are primarily activated according to Michiuchi et al. (2009).


**Table 5.1:** Parameters used for the crystal plasticity model, see Schäfer et al. (2019a).

We collect the three parameters we wish to determine into a single three-component vector

$$x = \begin{bmatrix} \tau\_c, A, B \end{bmatrix}^T \in \mathbb{R}^3,\tag{5.11}$$

which we would like to identify.

Following Herrera-Solaz et al. (2014) and Schäfer et al. (2019a), we consider the objective function

$$f^{h}(x) = \sqrt{\frac{1}{S} \sum\_{s=1}^{S} \left(\sigma\_{y,s}^{\text{exp}} - \sigma\_{y,s}^{\text{sim}}\right)^{2}},\tag{5.12}$$

quantifying the difference between a stress-strain hysteresis of low-cycle fatigue (LCF) experiments and the corresponding macroscopic stressstrain hysteresis of a micromechanical simulation. Here, denotes the number of time steps, exp refers to the experimentally measured stress in loading direction and sim stands for the homogenized stress in loading direction.

Following the experiments, we prescribe the macroscopic loading for the simulations to follow a triangular waveform in y-direction, see Fig. 5.2. In total, we simulate two cycles, where the each cycle lasts four seconds, and the second cycle serves as our simulation output. We use time increments of 0*.*01 seconds for computing the homogenized stress response, i.e., = 400 is used in this work.

Please note that we disregard strain-rate dependency of the considered material due to the cyclic stable behavior, see Schäfer et al. (2019a). In this work, we wish to minimize the mean of the error for different hysteresis, i.e., we consider

$$f(x) = \frac{1}{H} \sum\_{h=1}^{H} f^h(x). \tag{5.13}$$

Here, the term *ℎ* () corresponds to equation (5.12) and a fixed strain amplitude. In this article, we use the experimental data for the strain amplitudes = 0*.*35%, = 0*.*60% and = 0*.*90%, i.e., = 3. The corresponding hysteresis at half of the lifetime are shown in Fig. 5.3.

**Figure 5.3:** Polycrystalline stress-strain hysteresis of the martensitic 50CrMo4 steel at half of the lifetime used for the inverse optimization problem obtained by Schäfer et al. (2019a). (a) = 0*.*35%, (b) = 0*.*6%, (c) = 0*.*9%.

For a start, we study the necessary time-step size △. Indeed, exceedingly small time steps increase the computational cost, whereas too large time steps distort the computational results. We consider a very small time step of △ = 0*.*0001 s as our reference, and investigate the relative deviation of the error function (5.13) for time steps of increasing size, see Tab. 5.2. All simulations were carried out for the Armstrong-Frederick kinematic hardening model with the material parameters

found by Schäfer et al. (2019a). Based on these results, we select the time step △ = 0*.*01 s, as we consider the associated deviation of 1*.*26% in the error function negligible.


**Table 5.2:** Resulting error values for different time steps △.

For the parameter identification in a Bayesian optimization framework, diagrammatically represented in Fig. 5.1, we consider a two-stage process. First, we start with running BO on a large box in parameter space, to obtain a rough estimate for the parameters. Based on these results, a second optimization with tighter boundaries is carried out for fine-tuning the parameters. Both steps are described in further detail in the following subsections.

We use the GPy implementation (GPy, 2012) as a Gaussian-process regressor for the method described in Section 5.2. For the acquisition functions and the associated optimization we used an implementation by the Bosch Center for Artificial Intelligence (BCAI) based on Python 3.7, numpy (van der Walt et al., 2011) and scipy (Virtanen et al., 2020).

We use the maximum number of function evaluations as a stopping criterion within the BO framework. We set this value to 160, involving ten random initialization points and 150 BO steps, unless otherwise specified. For the evolutionary algorithms, included as a benchmark, we use the commercial software (Dynardo, 2019).

### **5.3.2 Parameter identification with a reasonably large search space**

For this setting, the search space is given by the parameter bounds listed in Tab. 5.3. We initialize Bayesian optimization with ten points from a Latin-Hypercube sampling. The "Matern5/2" kernel (5.5) is chosen as a kernel function. The lower confidence-bound acquisition-function (5.7) is employed with an exploration margin of = 2*.*0. Thus, over-exploitation and getting stuck in a local maximum are avoided. The rationale for choosing this specific acquisition function and exploration margin, is discussed in Section 5.3.3.


**Table 5.3:** Permissible parameter space for the large search space problem.

First, we investigate the parameter optimization for the AF model and investigate the smallest error, i.e., the smallest value of ,

$$f^\*: N \mapsto \min\_{i=1}^N f(x\_i),\tag{5.14}$$

versus function evaluations . Fig. 5.4a shows this quantity for a single BO run using a maximum number of function evaluations set to 400. Objective values above 80*.*0 MPa were cut off, as these are obtained during the first ten function evaluations in the Latin-Hypercube sampling.

Within this single run, BO reaches a minimum value of 48*.*80 MPa after 171 function evaluations. However, a comparable error of 49*.*9 MPa is reached after 81 evaluations.

BO reaches a strictly positive objective value. This indicates that the

material model described in Section 2.2.2 may only approximate the experimentally determined stress-strain hysteresis with a non-vanishing error. This effect may be systematic or result from stochastic fluctuations certainly present for the different experimentally determined hysteresis. The Bayesian optimization framework relies upon randomness both for the initialization and the optimization of the acquisition function. To evaluate the influence of this randomness, we restarted Bayesian optimization ten times. Fig. 5.4b shows the mean of ten BO runs in terms of the best error encountered versus number of function evaluations. Additionally, we computed the 95% confidence interval (95% CI) using an unbiased population formula for the standard deviation and Student's -distribution (Student, 1908). The bounds of the confidence interval are indicated by shaded areas centered at the mean value.

Please note that the number of evaluations used by BO is limited to 160, including those used for the initialization.

**Figure 5.4:** Best error versus function evaluations for the AF model using BO. (a) Single run, (b) Ten runs.

BO reaches a mean error of 47*.*83 MPa after the allowed 160 function evaluations with a 95% confidence interval of [46*.*40 MPa, 48*.*74 MPa]. A comparable mean error of 48*.*78 MPa and a 95% confidence interval [47*.*10*,* MPa, 48*.*76 MPa] is reached after 121 function evaluations, corresponding to 111 BO steps excluding evaluations used for initialization. Thus, we observe that repeating the BO process leads, after a sufficient number of iterations, to rather similar results as for a single run.

In addition to the AF model, we investigate the Ohno-Wang model in this work to demonstrate the general robustness and reliability of Bayesian optimization in this setting. The result for one BO run with a maximum number of 400 evaluations is shown in Fig. 5.5a.

Qualitatively, similar conclusions may be drawn as for the AF model. To address the issue of randomness we use ten BO runs with different initialization and evaluate the results, see Fig. 5.5b.

**Figure 5.5:** Best error versus function evaluations for the AF model using BO. (a) Single run, (b) Ten runs.

The observed behavior parallels the optimization for the AF model, where the mean error improves significantly within the first 100 function evaluations and stagnates afterwards. An important difference to the AF model are the final mean and the confidence bounds. For all considered optimization runs the mean best error is 43*.*94 MPa, which is lower than

the mean reached using AF. The 95% confidence interval is computed as [42*.*71 MPa, 45*.*17 MPa]. The broader 95% confidence interval for the OW model compared to the one obtained for the AF-model optimization indicates a higher dispersion of BO for the OW model, i.e., a higher influence of the initial sampling.

To gain insight into the shape of the energy landscape for the considered black-box function, we compute the difference between two consecutive minimum values of \* at the -th function evaluation

$$
\Delta f\_j^\* = f\_j^\* - f\_i^\* \quad \text{with} \quad i, j = 1 \ldots N,\tag{5.15}
$$

as well as the Euclidean distance between a parameter set (see equation (5.11)) and the previous vector −1. For proper dimensioning, we normalize each distance parameter-wise by the parameters \* , for which a minimum objective value was found, see Tab. 5.4, via

$$||x\_i - x\_{i-1}|| = \sqrt{\sum\_{d=1}^{D} \left(\frac{x\_{i,d} - x\_{i-1,d}}{x\_d^\*}\right)^2}.\tag{5.16}$$

Here, denotes the number of dimensions of the input vector (5.11), i.e., = 3 for this work. Let us denote by

$$\overline{\Delta f\_j^\*} = \frac{1}{R} \sum\_{r=1}^R \Delta f\_{j,r}^\* \tag{5.17}$$

$$\overline{d\_j} = \frac{1}{R} \sum\_{r=1}^{R} \left| \left| x\_j^r - x\_{j-1}^r \right| \right| \tag{5.18}$$

the averaged objective values and mean distances over total optimization runs at the -th function evaluation.


5 Identifying material parameters in crystal plasticity by Bayesian optimization

**Table 5.4:** Parameters used for normalization within the large search space setting.

First, we consider the Armstrong-Frederick kinematic hardening model, see Fig. 5.6a. The first ten function evaluations, corresponding to sampling of the given parameter space, allow the algorithm to find a promising starting point. Similar to Fig. 5.4b, for the first ten steps, the mean change in best error (encountered so far) is very high. After these initial steps, BO reduces the mean parameter distance while navigating towards a smaller error, i.e., the already gained knowledge about in the vicinity of an encountered (local) minimum is exploited and the uncertainty associated to this region is reduced in the process. Thereafter, it becomes more attractive to sample points with a higher degree of uncertainty, encoded by a higher variance. This exploration phase is indicated by the mean parameter distance oscillating around a fixed value at about 50 -112 function evaluations. The error is reduced frequently by small margins, up to about 112 function evaluations. After that, both the mean parameter distance and the mean difference in best error decrease further, as no further improvement of the objective value is made, see also Fig. 5.4b.

After about 120 function evaluations, corresponding to 110 BO steps, the exploitation tendency of the considered acquisition function is apparent. Despite decreasing the mean parameter distance,i.e., exploiting the already gained knowledge about the function, the best error is not decreased significantly. Indeed, the used acquisition function comes with an intrinsic trade-off between exploration and exploitation by sampling the point where the negative mean minus the variance (multiplied by ) is largest.

For the Ohno-Wang kinematic hardening model, the considered metrics

are shown in Fig. 5.6b. Sampling the parameter space reduces the encountered best error significantly, as we have seen for the Armstrong-Frederick model. Compared to the AF model, shown in 5.6a, is the mean relative parameter distance is smaller. This indicates a higher degree of smoothness in the objective function.

After the initialization phase, the error is reduced, exploiting the already gained knowledge of the objective function . The mean parameter distance decreases, i.e., neighboring points are being evaluated. After 50 function evaluations, the improvement of the best error gets less pronounced, resulting in the plateau apparent in Fig. 5.5b. After about 80 evaluations, a step with large mean parameter distance occurs, resulting in a subsequent decrease of the best error.

Interpreted in terms of the energy landscape of the objective , BO investigates a neighborhood with promising points first. After exploiting this area (up to about 64 iterations), the exploitation-exploration trade-off favors reducing the uncertainty. This results in a large mean parameter distance at about 80 function evaluations, allowing the algorithm to find a new area where the objective may be decreased further.

**Figure 5.6:** Comparison of the mean best error difference and parameter distance found by Bayesian optimization for the different kinematic hardening models. (a) Armstrong-Frederick, (b) Ohno-Wang.

For both single optimization runs shown in Fig. 5.4a and Fig. 5.5a, the resulting stress-strain hysteresis is shown in Fig. 5.7 besides their experimental counterparts. The results for, both, the Armstrong-Frederick and Ohno-Wang model, show good agreement with the experiments at medium and high strain amplitudes, i.e., = 0*.*6% and = 0*.*9%. These observations may be quantified when normalizing the root of the mean squared distance between experiment and simulation, i.e., equation (5.12), by to the root of the mean squared experimental values

$$\frac{f^h}{\sqrt{\frac{1}{N}\sum\_{i=1}^N \sigma\_{y,i}^{\text{exp}}}}.\tag{5.19}$$

For the Armstrong-Frederick model, the relative difference is 4*.*29% and 6*.*81% for = 0*.*6% and = 0*.*9%, respectively. We obtain similar values for the Ohno-Wang model, namely 3*.*56% and 6*.*48%. For the smallest strain amplitude considered, i.e., = 0*.*35%, the relative error is larger, with 22*.*60% and 20*.*41% for the Armstrong-Frederick and Ohno-Wang model, respectively. These errors may be reduced either by considering

a more complex microstructure or by working with more advanced models. For the work at hand, we focus on the optimization method used for the parameter identification.

**Figure 5.7:** A comparison of stress-strain hysteresis at different loading amplitudes with parameters identified by BO to experimental results. (a) = 0*.*35%, (b) = 0*.*6%, (c) = 0*.*9%.

#### **5.3.3 On the choice of acquisition function**

In this section, we elaborate on our choice of acquisition function and exploration margin. For the large search-space problem, we compare two different exploration margins for each of the considered acquisition functions and . For these four combinations, the mean best error of ten BO runs is shown in Fig. 5.8a. Except for with = 1*.*0, all combinations lead to a comparable minimum value.

**Figure 5.8:** Comparison of the mean best error using BO with different acquisition functions and exploration margins. (a) Entire range of evaluations, (b) Last 100 evaluations.

Additionally, we look at the maximum best error (i.e., the worst case) for all ten optimization runs versus the number of function evaluations, shown in Fig. 5.9a. As for the behavior of the mean best error, shown in Fig. 5.8 for with = 1*.*0, the maximum best error does not decrease to a comparable level. The other three considered cases behave similarly. To examine these cases in more detail, we investigate at the last 100 evaluations in Fig. 5.9b. BO using with an exploration margin of two reaches its minimum objective value of 49*.*45 MPa after 119 function evaluations for the considered case. The worst case for expected improvement and = 0*.*1, a stationary value of 51*.*36 MPa, is reached after 106 function evaluations. The reached value is not improved using more function evaluations. For an exploration margin of 0*.*05, a stationary value of 49*.*31 MPa is reached for the worst case considered. This is slightly smaller than the error obtained using and = 2*.*0, but requires 26 additional evaluations, which is why we chose with = 2*.*0 as our acquisition function and exploration margin.

**Figure 5.9:** Comparison of the worst case obtained best error using BO with different acquisition functions and exploration margins. (a) Entire range of evaluations, (b) Last 100 evaluations.

### **5.3.4 Fine-tuning material parameters with small parameter space**

In section 5.3.2, we were concerned with optimizing parameters when little is known about the optimum. We may work on a smaller parameter space if additional information is known, for instance based on expert knowledge or a previous optimization run for a similar problem. For determining the crystal plasticity parameters of the single crystal, we will use the bounds of Tab. 5.5, which we chose to be tighter than the previously considered bounds. Again, we choose our acquisition function to be with an exploration margin of = 2*.*0, as conclusions very similar to those in section 5.3.3 may be drawn the fine-tuning case, as well.


**Table 5.5:** Permissible parameter space for the fine-tuning problem.

The BO workflow is initialized by ten points obtained from a Latin-Hypercube sampling. Then, BO is run for a fixed set of 150 function evaluations, summing up to a total of 160. The results of ten optimization runs are shown in Fig. 5.10a and 5.10b for the Armstrong-Frederick and Ohno-Wang model, respectively.

**Figure 5.10:** Mean best error and corresponding 95% confidence interval versus function evaluations of the Bayesian optimization in the fine-tuning setting. (a) Armstrong-Frederick, (b) Ohno-Wang.

BO reaches a mean minimum error-value of 46*.*72 MPa for the AF model. The 95% confidence interval at end of optimization is given by [45*.*66 MPa, 47*.*78 MPa], which is tighter than for the results of the large search-space problem. After 80 evaluations, BO reaches a comparable mean error of 46*.*92 MPa with a 95% confidence interval of [45*.*76 MPa, 48*.*08 MPa].

For the optimization of the OW kinematic hardening model, Fig. 5.10b shows the mean and the 95% confidence interval of the smallest error versus function evaluations. The mean minimum-value of 41*.*69 MPa is reached after 136 evaluations and a comparable mean error of 41*.*695 MPa is reached after 73 evaluations. Furthermore, each of the optimization runs arrives at almost the same error, emphasized by the narrow 95% confidence interval s, [41*.*6851 MPa, 41*.*6854 MPa] and [41*.*52 MPa, 42*.*17 MPa] for 136 and 73 function evaluations, respectively.

As we did for the large search-space problem, we look at, both, the mean difference between two consecutive best encountered values of , and, the mean normed parameter distance (see equations (5.17) and (5.18)). As normalization, we use those parameters for which the minimum function value was found, see Tab. 5.6.


**Table 5.6:** Parameters used for normalization for the fine-tuning setting.

These metrics are shown for the AF model in Fig. 5.11a. BO systematically reduces the step size, with the main change made between 0 and 60 function evaluations. This corresponds to the observations made in Fig. 5.10a, where after about 60 function evaluations there is no substantial change in the mean best error encountered as well as for the 95% confidence interval . Fig. 5.11a also shows the trade-off between exploration and exploitation. This is indicated by peaks in parameter distance after about 80 function evaluations without decreasing . The mean distance between parameters is smaller within the more restrained boundaries of Tab. 5.5. This explains why a comparable or even better solution is achieved in shorter time than when optimizing with the

(larger) bounds given in Tab. 5.3.

For the fine-tuning of the OW model, see Fig. 5.11b, similar observations can be made. In contrast to optimizing the AF model, the mean change in best error, shown in Fig. 5.10b, is almost zero after 80 function evaluations. As for the large-space problem, the mean parameter distance is smaller than for the case of optimizing the Armstrong-Frederick model. However, for the case at hand, this difference is less pronounced.

**Figure 5.11:** Comparison of mean error and parameter distance obtained with BO for (a) the Armstrong-Frederick and (b) the Ohno-Wang kinematic hardening models within the fine-tuning setting.

#### **5.3.5 Comparison to the state-of-the-art**

In this section, we compare the proposed Bayesian approach to alternative optimization algorithms which demonstrated their power for similar tasks, in particular inversely identifying material parameters. In addition to the evolutionary approach used by Schäfer et al. (2019a), we investigate the Multilevel Coordinate Search (MCS) (Huyer and Neumaier, 1999) and Stable Noisy Optimization by Branch and FIT (SNOB-

FIT) (Huyer and Neumaier, 2008). The latter two derivative-free optimization methods turned out to be particularly powerful for lowdimensional problems in the expressive review article of Rios and Sahinidis (2013). The workflow used by Schäfer et al. (2019a) proceeds as follows. First, a Latin-Hypercube-Sampling, consisting of 200 points in the proposed parameter space, is performed. Subsequently, an evolutionary algorithm is run on the parameter space, initialized using the most promising results from the sampling step. MCS and SNOBFIT do not require an additional initialization, and we use the recommended default parameters (Huyer and Neumaier, 2008; 1999) for both algorithms.

We compare the smallest error encountered versus function evaluations for a single optimization run of the AF model in the large-search-space setting, see Fig. 5.12a, where we constrain the maximum number of function evaluations to 400. For a start, we observe that all algorithms reach a non-zero function value, in agreement to previous observations. The minimum error obtained by the evolutionary approach of 48*.*20 MPa is reached after 332 function evaluations. A comparable value of 49*.*15 MPa is encountered after 270 evaluations. Please note that the sampling provides a suitable initialization, and, for the first 100 function evaluations, there is a steady decrease in the objective value . MCS reaches a similar error of 48*.*69 MPa after 400 evaluations, with an error of 49*.*16 MPa reached after evaluating the cost function 237 times. Among all considered optimization algorithms, SNOBFIT reaches the lowest error value of 43*.*72 MPa after 201 function evaluations.

All considered algorithms, except for MCS, rely to some degree on a random initialization (Huyer and Neumaier, 2008; Brochu et al., 2010). Therefore, we wish to asses the influence of randomness on the overall optimization results. In Fig. 5.12b, the mean value of ten optimization runs is shown. Additionally, we visualize the computed 95% confidence interval via shaded areas, centered at the mean. For the SNOBFIT algorithm and these multiple runs, we set the maximum number of function evaluations to 160, as for the proposed Bayesian strategy. As

MCS is a deterministic algorithm, in all of the following we will only show the results of a single optimization run.

**Figure 5.12:** Comparison of the evolutionary approach, Bayesian optimization, SNOBFIT and MCS for the Armstrong-Frederick kinematic hardening model in the large search space setting. (a) Single run, (b) Ten runs.

After 100 evaluations, the evolutionary approach reaches a mean minimum function value of 56*.*31 MPa with a 95% confidence interval of [52*.*73 MPa, 59*.*90 MPa], and, after 160 evaluations arrives at a mean and 95% confidence interval of 53*.*42 MPa and [50*.*61 MPa, 56*.*23 MPa], respectively. The smallest mean value of 48*.*71 MPa is reached after evaluating the cost function 394 times, with a comparable value of 49*.*97 MPa reached after 376 evaluations. Using the SNOBFIT algorithm results in a smaller mean error of 52*.*66 MPa at 100 function evaluations, but a larger 95% confidence interval of [47*.*14 MPa, 58*.*19 MPa]. The mean best error after 160 function evaluations is lower for SNOBFIT than for BO and the evolutionary algorithm, namely 47*.*47 MPa while the 95% confidence interval of [44*.*15 MPa, 50*.*78 MPa], obtained by

using SNOBFIT, is larger than for the other two methods.

As our next step, we compare BO to the evolutionary approach, MCS and SNOBFIT for optimizing the unknown parameters of the OW kinematic hardening model in the large-space setting. In Fig. 5.13a, we plot the best error vs. iterations for each algorithm and a single optimization run. Bayesian optimization, the evolutionary approach and SNOBFIT reach similar error values of 44*.*21 MPa, 45*.*74 MPa and 43*.*72 MPa, respectively. Comparable error values are obtained by BO, SNOBFIT and the evolutionary framework after evaluating the cost function about 120, 290 and 340 times, respectively. Using Multilevel Coordinate Search leads to a minimum best error of 48*.*69 MPa, which is reached after 392 function evaluations.

To investigate the influence of the random initialization on optimizing the OW model, we use the mean of ten optimization runs and visualize the results, with the corresponding 95% confidence interval as shaded areas, in Fig. 5.13b for BO, SNOBFIT and the evolutionary approach.

Combining Latin-Hypercube sampling and the evolutionary algorithm produces, after 400 function evaluations, a mean error of 51*.*10 MPa, and a 95% confidence interval of [47*.*76 MPa, 54*.*40 MPa]. After evaluating the function 160 times, the mean best error using the evolutionary approach is 55*.*30 MPa with 95% confidence interval [53*.*19 MPa, 57*.*40 MPa].

For SNOBFIT, the mean best error of 49*.*77 MPa after 160 function evaluations lies between the results of BO and the evolutionary approach, whereas the 95% confidence interval of [44*.*73 MPa, 54*.*81 MPa] is larger than for the two other methods.

**Figure 5.13:** Comparison of the evolutionary approach, Bayesian optimization, SNOBFIT and MCS for the Armstrong-Frederick kinematic hardening model in the large search space setting. (a) Single run, (b) Ten runs.

In Fig. 5.14, we compare the considered algorithms for the refined bounds given in Tab. 5.5, for both kinematic hardening models.

For optimizing the parameters of the AF model, SNOBFIT provides the lowest mean error of 44*.*16 MPa with [42*.*99 MPa, 45*.*34 MPa] as a 95% confidence interval after 160 function evaluations. We observe comparable values after letting SNOBFIT evaluate the function 98 times, e.g., 44*.*53 MPa for the mean error with a 95% confidence interval of [43*.*32 MPa, 45*.*74 MPa]. Using an evolutionary algorithm to optimize the AF parameters leads to a mean minimum error of 47*.*29 MPa, which is close to the mean value reached by BO. The 95% confidence interval at 160 function calls of the evolutionary algorithm computes to [45*.*47 MPa and 49*.*10 MPa], which is larger than for BO or SNOBFIT. To reach similar mean error values of 47*.*50 MPa and 46*.*85 MPa, the evolutionary algorithm and BO need 153 and 87 function evaluations, respectively. Among the considered algorithms, MCS reaches the lowest best error,

already in the first 13 iterations with 50*.*807 MPa. This error level is reached within the first function evaluations. Yet, MCS falls short of reducing the error significantly for subsequent function evaluations. Indeed, after evaluating the function 160 times, the error is 48*.*70 MPa. Next, we compare the mean best error encountered when optimizing the parameters for the Ohno-Wang kinematic hardening model using the evolutionary approach, BO or SNOBFIT, see Fig. 5.14b, together with best error over iterations for the MCS method. Similar conclusions for the behavior of SNOBFIT may be drawn as for the BO approach, see Sec. 5.3.4. The mean best error and 95% confidence interval after 160 function evaluations are 43*.*10 MPa and [42*.*60 MPa, 43*.*60 MPa], respectively. Using the evolutionary algorithm produces a mean minimum error of 49*.*31 MPa after 133 function evaluations with a 95% confidence interval of [48*.*78 MPa, 49*.*84 MPa]. Almost the same mean best error of 50*.*71 MPa with a 95% confidence interval of [48*.*79 MPa, 52*.*64 MPa] is reached after 13 function evaluations and not improving up to 109 evaluations. For MCS, the error is, similar to optimizing the AF model, the lowest in the first 10 optimization steps. After evaluating the cost function 160 times, the error reaches a similar value as reached by BO with 42*.*29 MPa.

**Figure 5.14:** Comparison of the evolutionary approach, Bayesian optimization, SNOBFIT and MCS for both kinematic hardening model in the fine-tuning setting. (a) Armstrong-Frederick model, (b) Ohno-Wang model.

Let us summarize the findings of this section for each algorithm individually. The evolutionary algorithm consistently required the highest number of function evaluation among the considered algorithms to reach a specific error level. For the OW model and both the large and the small search space, the evolutionary algorithm lacked significantly behind the other algorithms. With the exception of the large search space and the AF model, MCS provided a similar error value as BO and the evolutionary algorithm. MCS showed inferior performance for the OW model and the large search space. In the case of fine tuning, MCS turned out to be worst for the Armstrong-Frederick model, but second best for Ohno-Wang. Still, there are practical benefits of the MCS method. First and foremost, it is a *deterministic* algorithm, dispensing with the influence of randomness. Secondly, in the fine-tuning setting MCS provided very good results already for the first few iterations. Thus, when the number of evaluations is small and prescribed beforehand, MCS is the method of choice.

SNOBFIT is good for the fine-tuning case, in particular for the AF model. In this case, SNOBFIT produces a similar dispersion as BO. However, for the large search space, the values produced by SNOBFIT are characterized by a rather large dispersion. Directly comparing SNOBFIT and BO, the former consistently requires more function evaluations to reach the same error level.

The proposed Bayesian optimization framework produces a mean which is comparable to SNOBFIT in the setting of the large search space. However, BO is both consistently faster and characterized by a smaller dispersion. Only for fine-tuning AF, BO is a little worse than SNOBFIT in terms of the reached best error. For fine-tuning the more complex Ohno-Wang model, BO is fastest and leads to the best error.

To conclude, the presented BO framework provides a strongly competitive solution, independent of the problem. Thus, it turns out to be the algorithm of choice for industrial applications, in particular parameter identification for an unknown material.

#### **5.3.6 Industrial-scale polycrystalline microstructure**

Previously, we determined the optimum parameters for a simplified volume element, see Fig. 5.2, which was selected with extensive studies of the proposed Bayesian optimization approach in mind. In this section, we investigate a more complex microstructure to demonstrate the capabilities of the approach. For this we use a synthetic microstructure, generated by the method discussed in Kuhn et al. (2020), consisting of 100 grains (with equal volume). Furthermore, each grain is randomly assigned a crystalline orientation, s.t. the overall orientation is isotropic, see Fig. 5.16a, where the coloring indicates the corresponding orientation. We furnish this volume element with periodic boundary conditions, discretize it by 64<sup>3</sup> voxels and use the FFT-based solver FeelMath (Fraunhofer ITWM, 2021; Wicht et al., 2020a;b) to speed up

computations (Rovinelli et al., 2020). Both, the Armstrong-Frederick and the Ohno-Wang kinematic hardening model will be considered.

Following Sedighiani et al. (2020), we define the bounds of the feasible parameter space in terms of the previously found optima, with a range of 1*.*5 × \* and 0*.*5 × \* , see Tab. 5.7 and Tab. 5.8. We retain our choice of acquisition function and exploration margin. To keep the computational cost reasonable, we restrict the maximum number of iterations to 50.


**Table 5.7:** Parameter space for the AF model following Sedighiani et al. (2020).


**Table 5.8:** Parameter space for the OW model following Sedighiani et al. (2020).

The best error versus number of evaluations is shown in Fig. 5.15a for the Armstrong-Fredrick model. Guided by a Latin-Hypercube-Sampling, the first ten function evaluations result in an error of 44*.*10 MPa. This error is further reduced to 43*.*68 MPa after 20 additional function evaluations. The minimum error of 41*.*20 MPa is reached after evaluating the function 33 times. The obtained error is smaller than in Section 5.3.4. This may be a result of the increased complexity of the microstructure, which is closer to the experimental setup encompassing a large number of grains. Furthermore, the obtained parameter values, shown in Tab. 5.9 lie outside of the bounds considered in the previous sections, see Tab. 5.3 and Tab. 5.5.

**Figure 5.15:** Smallest error versus iterations for the optimization of kinematic hardening parameters and critical resolved shear stress using the microstructure with increased complexity. (a) Armstrong-Frederick model, (b) Ohno-Wang model.

Similar conclusions may be drawn for the Ohno-Wang kinematic hardening model, where the smallest error versus number of function evaluations is shown in Fig 5.15b. After ten function evaluations following a Latin-Hypercube-Sampling, the smallest error is 62*.*05 MPa. After taking one BO step, this error is reduced to 41*.*35 MPa. The final minimum error of 38*.*21 MPa is reached after 27 function evaluations, i.e., after 17 steps proposed by Bayesian optimization.


**Table 5.9:** Resulting parameter sets for the optimization with a more complex microstructure model.

To ensure that the smaller error for the complex microstructure is not a consequence of the different parameter bounds, we return to the simplified microstructure, see Fig. 5.2, and run Bayesian optimization on the constrained parameter spaces, see Tab. 5.7 and Tab. 5.8. For the Armstrong-Frederick model, BO reaches a mean error of 44*.*69 MPa an a [44*.*69 MPa, 44*.*70 MPa] 95% confidence interval after 50 function evaluations. After the same number of evaluations, the resulting mean error and 95% confidence interval are 41*.*78 MPa and [41*.*59 MPa, 41*.*97 MPa], respectively. Thus, even for the constrained bounds, the error corresponding to the more complex microstructure shown in Fig. 5.16a improves upon the error obtained for the simplified structure.

Concerning the computational cost, the simulation with a time step of △ = 0*.*01s and 16 CPUs takes roughly 40 minutes for the largest strain amplitude of = 0*.*9%. Thus, for the chosen maximum number of 50 function evaluations, the total Bayesian parameter identification takes 33 hours. The minimum error of 38*.*21 MPa is found after a computation time of 18 hours in total for the Ohno-Wang kinematic hardening model. As for the Armstrong-Frederick model a total computation time of 22 hours is needed to find a solution with error 41*.*20 MPa.

Using a more complex synthetic microstructure, the error is further reduced, at the expense of an increased computational effort. Still, the Bayesian-optimization approach limits the required total time to less than a day. Using a more complex representation of the microstructure offers also permits us to gain deeper insights into the deformation behavior of the microstructure, for instance in terms of the accumulated plastic slip, a mesoscopic indicator for the plastic deformation. Furthermore, the obtained optimum parameter set may enter further microstructure simulations.

For = 0*.*9% and = 8*.*0s, the distribution of accumulated plastic slip, see Wulfinghoff et al. (2013),

$$\dot{\gamma} = \sum\_{\eta=1}^{N\_S} \left| \dot{\gamma}\_{\eta} \right|, \tag{5.20}$$

is shown in Fig. 5.16b and Fig. 5.16c for the Armstrong-Frederick and Ohno-Wang kinematic hardening model, respectively. In both cases, the accumulated plastic slip is concentrated in specific grains.

**Figure 5.16:** Resulting distribution of the accumulated plastic ˙ slip for = 0*.*9% at the end of simulation. (a) Morphology, (b) AF model, (c) OW model.

As we used the same geometric microstructure for each of the hardening models, we observe that, independent of the employed kinematic hardening model, the orientation of each grain plays a crucial role in determining the deformation behavior. Apparently, grains with higher accumulated plastic slip ˙ are oriented in such a way that the applied macroscopic load has maximum effect, thus resulting in a higher amount of plastic deformation.

Moreover, the distribution of accumulated plastic slip provides some insight into the advantages of the OW model over the AF model in producing smaller errors which is not apparent from the resulting hysteresis shown in Fig. 5.17. For the Ohno-Wang kinematic hardening model, there are more regions with high accumulated plastic slip ˙ in Fig. 5.16c than for the Armstrong-Frederick model in Fig. 5.16b. This more pronounced plasticity seems to match the experimental behavior better than the results computed using the Armstrong-Frederick kinematic hardening model. The relative error defined in equation (5.19) is

considerably lower for the small strain amplitude = 0*.*35% for both kinematic hardening models with 11*.*65% and 8*.*65% for the Armstrong-Frederick and Ohno-Wang model, respectively. The relative errors for the medium strain amplitude are increased to 11*.*20% for the Armstrong-Frederick and 11*.*08% for the Ohno-Wang model, whereas the error for = 0*.*9% is also reduced to 4*.*59% and 5*.*16% for the respective models. As the small strain amplitude = 0*.*35% is most interesting for fatigue applications, the benefits of using the more complex structure, see Fig. 5.16a, sinstead of the simplified microstructure, see Fig 5.2, become apparent.

**Figure 5.17:** Comparison of stress-strain hysteresis determined by BO to experimental results at different loading amplitudes. (a) = 0*.*35%, (b) = 0*.*6%, (c) = 0*.*9%.

## **5.4 Conclusion**

This work was dedicated to identifying constitutive parameters for a small strain crystal plasticity constitutive law incorporating a kinematic hardening model to capture the behavior under cyclic loading conditions, a topic which received attention recently from, both, a scientific and applied perspective (Sedighiani et al., 2020; Shahmardani et al., 2020; Shenoy et al., 2008; Prithivirajan and Sangid, 2018; Guery et al., 2016; Hochhalter et al., 2020; Kapoor et al., 2021; Pagan et al., 2017; Bandyopadhyay et al., 2020; 2019; Rovinelli et al., 2018).

We investigated a Bayesian optimization approach for the rapid and flexible calibration of the single-crystal constitutive parameters entering micromechanical simulations. The methodology is data-driven and compares the predicted stress-strain hysteresis to given polycrystalline experimental data. More precisely, we proposed using Bayesian optimization with Gaussian processes to build up a surrogate model of the optimization problem at hand and employed the upper confidence bound using an exploration margin of two for selecting the next iterate. Our numerical investigations demonstrated that the proposed optimization framework reliably and efficiently determines suitable material parameters for our purposes. For instance, these moduli may enter a simulation predicting the initiation of fatigue cracks, see Schäfer et al. (2019a). We investigated both a large and a comparatively small space of admissible parameters, simulating whether prior knowledge on the optimum is available or not. The proposed framework is able to handle the challenges of the large search space, and is sped up for the smaller admissible set.

By investigating both the changes in parameter vector and in the objective value, we could gain an intuitive understanding for the selection strategy of Bayesian optimization, providing a suitable tradeoff between exploitation and exploration. Furthermore, we could get an insight into the general energy landscape and the effects of using different acquisition functions.

We compared the proposed Bayesian optimization framework to a representative Schäfer et al. (2019a) of the powerful class of genetic algorithms (Kapoor et al., 2021; Bandyopadhyay et al., 2020; Rovinelli et al., 2018) as well as two derivative-free optimization schemes recommended in the review article by Rios and Sahinidis (2013) for low-dimensional problems. The computational findings of section 5.3.5 demonstrated the

capabilities of the Bayesian optimization framework. Compared to the investigated evolutionary algorithm, BO turned out to be consistently faster, and featured a smaller dispersion. BO outperforms MCS mainly in the case of the large search space, whereas, for fine tuning, BO produces a smaller minimum achieved error for all cases, but requires more function evaluations to reach a comparable error level. Compared to SNOBFIT, BO produces a considerably smaller dispersion when considering a large search-space. BO requires fewer function evaluations, and provides similar or better results in case of fine-tuning.

We observed that, for fitting computed stress-strain hysteresis to their experimental counterparts accurately, using more complex microstructures is necessary. In particular, the hysteresis at a low strain amplitude, the one most relevant to fatigue applications, using a volume element with 100 complex grains led to more accurate results than the simplified volume element with 64 cubic grains.

To further improve the fitting quality, it appears wise to investigate whether increasing the complexity of the underlying material model or tweaking the parameters used for microstructure generation has a higher effects. Especially in the context of industrial applications and the digital twin paradigm (Tuegel et al., 2011), accurate computational representations of components and their microstructure are needed. Further research might be invested into optimizing the microstructure descriptors, which are used as input for synthetic microstructure generation (Quey and Renversade, 2018; Bourne et al., 2020; Kuhn et al., 2020). In particular, it may be possible to determine a set of minimal, thus most efficient, characteristics a microstructure has to possess in order to produce the desired macroscopic material behavior. Complementing research dedicated to microstructure representations, the proposed optimization scheme could be improved. Despite the power and flexibility we experienced in the applied context, decreasing the number of necessary simulation runs (via fewer function calls) would be very much appreciated, in particular combined with a clever

stopping criterion. Furthermore, research could be devoted to dedicated transfer strategies, i.e., building upon previous experience with similar microstructures, material models or experimental data (Golovin et al., 2017). Last but not least, it appears worthwhile to extend the Bayesian optimization framework to uncertainty quantification, see Bandyopadhyay et al. (2019).

### **Chapter 6**

# **Summary and conclusions**

To ensure the safety of industrial components made from polycrystalline metals, it is imperative to assess their reliability with respect to cyclic loading conditions, as fatigue is still one of the major root cause of failure. Heterogeneities on the microscale are sources of fatigue crack initiation. Thus, the material's microstructure strongly impacts the fatigue behavior of a component. Typically, these influences are investigated by means of experiments. However, fatigue experiments are time and cost intensive. Simulating the mechanical behavior on a microscopic scale allows predicting fatigue lifetimes, taking the microstructures influence into account and reducing experimental effort.

In this work we adressed two topics required for micromechanical simulations of polycrystalline metals: Providing representative volume elements (RVEs) and identifying parameters of the crystal plasticity model in use. With the goal of reducing the overall runtime of the micromechanical simulations, we proposed two methods to create synthetic polycrystalline microstructures which reduce the required RVE size. We used a Bayesian optimization approach to inversely identify the model parameters based on polycrystalline experiments, reducing the required simulations by about a factor of two.

We considered Laguerre tessellations as models for polycrystalline microstructures in Chap. 3. Based on the formulation as an convex optimization problem, Bourne et al. (2020) proposed a method to quickly compute the weights of Laguerre tessellations, so that the computed cells have

prescribed volume fractions. We showed that using modern gradienttype solvers enables computing the weights within a few seconds. For all numerical examples considered, the non-monotone Barzilai-Borwein scheme provided the lowest run-times and therefore may serve as a general-purpose method. For the sake of shape regularity, it might be beneficial that the cells are *centroidal*, i.e., their seeds coincide with their centroids. We applied Anderson acceleration to Lloyd's algorithm used by Bourne et al. (2020), which considerably speeds up the computations. With the algorithm at hand we were able to compute shape regular Laguerre tessellations which reproduced an experimentally observed grain size distribution. Additionally, computational cells consisting of multiple phases can be generated, e.g., accounting for different phase compositions due to a heat treatment of the metal. Extending the method to other features of the microstructure, e.g., sphericity of the grains, might be an interesting step to increase the physical realism of the generated synthetic microstructures.

As a next step, in Chap. 4, we developed a method for prescribing crystallographic orientations to individual grains of the polycrystalline microstructures. We used the coefficients of a Fourier series expansion of the crystallite orientation distribution function, to formulate an iterative update procedure for the crystallographic orientations in the volume element. The cost function was formulated as the difference between prescribed and realized so-called texture coefficients. We compared the proposed method to other state-of-the-art approaches in terms of required volume element size to reproduce prescribed effective behavior. We considered, both, a uniform as well as a textured orientation distribution and investigated linear as well as non-linear behavior. For a unique grain size distribution, we observed the proposed optimization method to require a smaller number of grains in the volume element, thus, reducing the size of the RVE. Extending our investigations to the case of a log-normal grain size distribution, we found the proposed "*T*exture coefficient *O*ptimization for *P*rescribing orientations" (TOP) method shows little sensitivity to the underlying grain sizes, in contrast to the other considered methods. To further improve the proposed method, it would be interesting if faster optimization schemes, e.g., a Newton method, could be applied.

With a framework, enabling the creation of high-fidelity microstructures at hand, we proposed a Bayesian optimization approach to determine suitable material model parameters in Chap. 5. More precisely, using an inverse optimization approach and Gaussian processes as surrogate models, we were able to determine *single* crystal plasticity parameter from *polycrystalline* experiments quickly and efficiently. Compared to the state-of-the-art and powerful derivative-free approaches, the Bayesian optimization proved to be the most robust and the fastest. Applying transfer-learning approaches could further improve the efficiency of the optimization, especially when facing new, but similar, materials for which the parameters are to be determined. Additionally, it might be interesting to investigate whether the approach can be used to determine *single* phase parameters using experimental results from a *multiphase* specimen, e.g., when it is not possible to manufacture specimens with a single phase.

Taking all the described results into account, this thesis supplied novel methods to create microstructural representations of polycrystalline materials and provided a fast and reliable optimization method to parameterize the underlying material model. The promoted methods in this thesis are expected to accelerate the development of digital twins to drive the digitalization of the manufacturing process of metallic components. In order to predict the lifetime of components under various loading conditions using a digital twin, suitable methods for predicting the mechanical behavior have to be used. The methods developed in this thesis enable creating micromechanical models, which might be used either in FE<sup>2</sup>/FE-FFT approaches (Spahn et al., 2014; Kochmann et al., 2016) or modern, Machine Learning based methods (Gajek et al., 2021).

# **Bibliography**

Adams, B., Olson, T., 1998. The mesostructure - property linkage in polycrystals. Progress in Materials Science 43 (1), 1–87.

Adams, B. L., Boehler, J. P., Guidi, M., Onat, E. T., 1992. Group theory and representation of microstructure and mechanical behavior of polycrystals. Journal of the Mechanics and Physics of Solids 40 (4), 723–737.

Allen, M. P., Tildesley, D. J., 1987. Computer Simulation of Liquids. Clarendon Press, Oxford.

Alpers, A., Brieden, A., Gritzmann, P., Lyckegaard, A., Poulsen, H. F., 2015. Generalized balanced power diagrams for 3D representations of polycrystals. Philosophical Magazine 95 (9), 1016–1028.

Anderson, D. G., 1965. Iterative procedures for nonlinear integral equations. Journal of the ACM 12 (4), 547–560.

Anderson, E., Bai, Z., Bischof, C., Blackford, S., Demmel, J., Dongarra, J., Du Croz, J., Greenbaum, A., Hammarling, S., McKenney, A., Sorensen, D., 1999. LAPACK Users' Guide, 3rd Edition. Society for Industrial and Applied Mathematics, Philadelphia.

Anderson, M., Grest, G., Srolovitz, D., 1989. Computer simulation of normal grain growth in three dimensions. Philos Magaz B 59 (3), 293–329.

Arnaudov, N., Kolyshkin, A., Weihe, S., 2020. Micromechanical modeling of fatigue crack initiation in hydrogen atmosphere. Mechanics of Materials 149, 103557.

Arts, R. J., 1993. A study of general anisotropic elasticity in rocks by wave propagation: Theoretical and experimental aspects. Ph.D. thesis, Institut français du pétrole.

Asaro, R. J., Needleman, A., 1985. Texture development and strain hardening in rate dependent polycrystals. Acta Metallurgica 33 (6), 923–953.

Asaro, R. J., Rice, J. R., 1977. Strain localization in ductile single crystals. Journal of the Mechanics and Physics of Solids 25 (5), 309–338.

Asmussen, S., Glynn, P. W., 2007. Stochastic Simulation: Algorithms and Analysis. Springer, New York.

Aurenhammer, F., 1991. Voronoi diagrams - a survey of a fundamental geometric data structure. ACM Computing Surveys 23, 345–405.

Aurenhammer, F., Hoffmann, F., Aronov, B., 1998. Minkowski-Type Theorems and Least-Squares Clustering. Algorithmica 20, 61–76.

Ball, J. M., Carstensen, C., 2015. Geometry of polycrystals and microstructure. arXiv e-prints 1507.05197, 1–7.

Bandyopadhyay, R., Prithivirajan, V., Peralta, A. D., Sangid, M. D., 2020. Microstructure-sensitive critical plastic strain energy density criterion for fatigue life prediction across various loading regimes. Proceedings of the Royal Society A 476 (2236), 20190766.

Bandyopadhyay, R., Prithivirajan, V., Sangid, M. D., 2019. Uncertainty Quantification in the Mechanical Response of Crystal Plasticity Simulations. JOM 71 (8), 2612–2624.

Bansal, R. K., Kubis, A., Hull, R., Fitz-Gerald, J., 2006. High-resolution three-dimensional reconstruction: a combined scanning electron microscope and focused ion-beam approach. Journal of Vacuum Science & Technology B 24 (2), 554 – 561.

Barbe, F., Decker, L., Jeulin, D., Cailletaud, G., 2001. Intergranular and intragranular behavior of polycrystalline aggregates. Part 1: FE model. International Journal of Plasticity 17, 513–536.

Bargmann, S., Klusemann, B., Markmann, J., Schnabel, J. E., Schneider, K., Soyarslan, C., Wilmers, J., 2018. Generation of 3D representative volume elements for heterogeneous materials: A review. Progress in Materials Science 96, 322–384.

Bari, S., Hassan, T., 2000. Anatomy of coupled constitutive models for ratcheting simulation. International Journal of Plasticity 16 (3-4), 381– 409.

Barzilai, J., Borwein, J. M., 1988. Two-point step size gradient methods. IMA Journal of Numerical Analysis 8, 141–148.

Becker, R., 1991. Analysis of texture evolution in channel die compression. 1. Effects of grain interaction. Acta Metallurgica et Materialia 39, 1211–1230.

Bertram, A., 2012. Elasticity and Plasticity of Large Deformations. Springer, Berlin Heidelberg, Heidelberg.

Bertram, A., Böhlke, T., Gaffke, N., Heiligers, B., Offinger, R., 2000. On the generation of discrete isotropic orientation distributions for linear elastic cubic crystals. Journal of elasticity and the physical science of solids 58 (3), 233–248.

Bishop, C. M., 2006. Pattern Recognition and Machine Learning. Springer, New York.

Bishop, J. F. W., 1953. VI. A theoretical examination of the plastic deformation of crystals by glide. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 44 (348), 51–64.

Biswas, A., Prasad, M. R. G., Vajragupta, N., Kostka, A., Niendorf, T., Hartmaier, A., 2020. Effect of Grain Statistics on Micromechanical Modeling: The Example of Additively Manufactured Materials Examined by Electron Backscatter Diffraction. Advanced Engineering Materials 22 (5), 1901416.

Blum, C., Merkle, D., 2008. Swarm Intelligence: Introduction and Applications. Springer, Berlin.

Böhlke, T., 2005. Application of the maximum entropy method in texture analysis. Computational Materials Science 32 (3-4), 276–283.

Böhlke, T., 2006. Texture simulation based on tensorial Fourier coefficients. Computers & Structures 84 (17-18), 1086–1094.

Böhlke, T., Bertram, A., 2001. Isotropic orientation distributions of cubic crystals. Journal of the Mechanics and Physics of Solids 49 (11), 2459– 2470.

Böhlke, T., Bertram, A., 2003. Crystallographic texture induced anisotropy in copper: An approach based on a tensorial Fourier expansion of the codf. In: Journal de Physique IV (Proceedings). Vol. 105.

Böhlke, T., Jöchen, K., Kraft, O., Löhe, D., Schulze, V., 2010. Elastic properties of polycrystalline microcomponents. Mechanics of Materials 42 (1), 11–23.

Böhlke, T., Lobos, M., 2014. Representation of Hashin–Shtrikman bounds of cubic crystal aggregates in terms of texture coefficients with application in materials design. Acta Materialia 67, 324–334.

Borbely, A., Biermann, H., Hartmann, O., 2001. FE investigation of the effect of particle distribution on the uniaxial stress–strain behaviour of particulate reinforced metal-matrix composites. Materials Science and Engineering: A 313 (1-2), 34–45.

Bourne, D. P., Kok, P. J. J., Roper, S. M., Spanjer, W. D. T., 2020. Laguerre tessellations and polycrystalline microstructures: A fast algorithm for generating grains of given volumes. Philosophical Magazine 100 (21), 2677–2707.

Bourne, D. P., Roper, S. M., 2015. Centroidal Power diagrams, Lloyd's algorithm and applications to optimal location problems. SIAM Journal on Numerical Analysis 53 (6), 2545–2569.

Boyd, S., Vandenberghe, L., 2004. Convex Optimization. Cambridge University Press, Cambridge.

Bravais, A., 1850. Les Systemes Formes par des Points Distribues Regulierement sur un Plan ou Dans Lespace. Journal de l'Ecole Polytechnique, XIX 1, 1–128.

Brenner, R., 2009. Numerical computation of the response of piezoelectric composites using Fourier transform. Physical Review B 79 (18), 184106.

Brieden, A., Gritzmann, P., 2012. On optimal weighted balanced clusterings: Gravity bodies and power diagrams. SIAM Journal on Discrete Mathematics 26 (2), 415–434.

Brisard, S., Dormieux, L., 2010. FFT-based methods for the mechanics of composites: A general variational framework. Computational Materials Science 49 (3), 663–671.

Brochu, E., Cora, V. M., De Freitas, N., 2010. A Tutorial on Bayesian Optimization of Expensive Cost Functions, with Application to Active User Modeling and Hierarchical Reinforcement Learning. arXiv preprint arXiv:1012.2599 .

Broyden, C. G., 1970. The convergence of a class of double rank minimization algorithms: 2. The new algorithm. Journal of Mathematical Analysis and Applications 6, 222–231.

Bunge, H.-J., 1982. Texture Analysis in Materials Science: Mathematical Methods. Cuvillier Verlag.

Burdakov, O., Dai, Y.-H., Huang, N., 2019. Stabilized Barzilai-Borwein method. Journal of Computational Mathematics 37 (6), 916–936.

Cailletaud, G., 1992. A micromechanical approach to inelastic behaviour of metals. International Journal of Plasticity 8 (1), 55–73.

Cavallini, F., 1999. The best isotropic approximation of an anisotropic Hooke's law. Bollettino di Geofisica Teorica ed Applicata 40 (1), 1–18.

Chakraborty, A., Eisenlohr, P., 2017. Evaluation of an inverse methodology for estimating constitutive parameters in face-centered cubic materials from single crystal indentations. European Journal of Mechanics-A/Solids 66, 114–124.

Chen, Y., Vasiukov, D., Gélébart, L., Park, C. H., 2019. A FFT solver for variational phase-field modeling of brittle fracture. Computer Methods in Applied Mechanics and Engineering 349, 167–190.

Chunlei, X., Nakamachi, E., Xianghuai, D., 2000. Study of texture effect on strain localization of BCC steel sheets. Acta Mechanica Solida Sinica 13 (2), 95–104.

Coleman, B. D., Noll, W., 1974. The Thermodynamics of Elastic Materials with Heat Conduction and Viscosity. In: The Foundations of Mechanics and Thermodynamics. Springer, Berlin, Heidelberg, Heidelberg, pp. 145–156.

Cottrell, A. H., 1954. Dislocations and Plastic Flow in Crystals. American journal of physics 22 (4), 242–243.

Courtney, T. H., 2005. Mechanical behavior of materials. Waveland Press, Long Grove.

Cox, D. D., John, S., 1992. A statistical method for global optimization. In: 1992 IEEE International Conference on Systems, Man, and Cybernetics. IEEE.

Cruzado, A., Gan, B., Jiménez, M., Barba, D., Ostolaza, K., Linaza, A., Molina-Aldareguia, J. M., Llorca, J., Segurado, J., 2015. Multiscale modeling of the mechanical behavior of IN718 superalloy based on micropillar compression and computational homogenization. Acta Materialia 98, 242–253.

Cruzado, A., Llorca, J., Escudero, J. S., 2020. Computational Micromechanics Modeling of Polycrystalline Superalloys: Application to Inconel 718. In: Ghosh, S., Woodward, C., Przybyla, C. (Eds.), Integrated Computational Materials Engineering (ICME). Springer, New York, pp. 127–163.

Dai, Y.-H., 2012. A perfect example for the BFGS method. Mathematical Programming 138 (1–2), 501–530.

Dai, Y. H., Liao, L. Z., 2002. R-linear convergence of the Barzilai and Borwein gradient method. IMA Journal of Numerical Analysis 22, 1–10.

Dassault Systèmes, 2018. ABAQUS. Version: 2018.

URL https://www.3ds.com/de/produkte-und-services/ simulia/produkte/abaqus/

Deka, D., Joseph, D. S., Ghosh, S., Mills, M. J., 2006. Crystal plasticity modeling of deformation and creep in polycrystalline Ti-6242. Metallurgical and Materials Transactions A 37 (5), 1371–1388.

Dembo, R. S., Eisenstat, S. C., Steihaug, T., 1982. Inexact Newton Methods. SIAM Journal on Numerical Analysis 19, 400–408.

Döbrich, C. M., Rau, C., Krill III, C. E., 2004. Quantitative characterization of the three-dimensional microstructure of polycrystalline Al-Sn using X-ray microtomography. Metallurgical and Materials Transactions A 35 (7), 1953 – 1961.

Dyck, A., Böhlke, T., 2020. A micro-mechanically motivated phenomenological yield function for cubic crystal aggregates. ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik 100 (4), e202000061.

Dynardo, 2019. optiSLang. Version: 3.3.5.

URL https://www.dynardo.de/software/optislang.html

Eisenlohr, P., Roters, F., 2008. Selecting a set of discrete orientations for accurate texture reconstruction. Computational Materials Science 42 (4), 670–678.

Engels, J. K., Vajragupta, N., Hartmaier, A., 2019. Parameterization of a non-local crystal plasticity model for tempered lath martensite using nanoindentation and inverse method. Frontiers in Materials 6, 247.

Ernesti, F., Schneider, M., Böhlke, T., 2020. Fast implicit solvers for phasefield fracture problems on heterogeneous microstructures. Computer Methods in Applied Mechanics and Engineering 363, 112793.

Eshelby, J. D., 1957. The determination of the elastic field of an ellipsoidal inclusion, and related problems. Proceedings of the royal society of London. Series A. Mathematical and physical sciences 241 (1226), 376– 396.

Evans, C., Pollock, S., Rebholz, L. G., Xiao, M., 2020. A proof that Anderson acceleration improves the convergence rate in linearly converging fixed point methods (but not in those converging quadratically). SIAM Journal on Numerical Analysis 58 (1), 788–810.

Fang, H.-R., Saad, Y., 2009. Two classes of multisecant methods for nonlinear acceleration. Numerical Linear Algebra with Applications 16, 197–221.

Farooq, H., Cailletaud, G., Forest, S., Ryckelynck, D., 2020. Crystal plasticity modeling of the cyclic behavior of polycrystalline aggregates under non-symmetric uniaxial loading: Global and local analyses. International Journal of Plasticity 126, 102619.

Fedorov, F. I., 1968. Theory of Elastic Waves in Crystals. Springer Science & Business Media.

Fernández, M. L., Böhlke, T., 2019. Representation of Hashin–Shtrikman Bounds in Terms of Texture Coefficients for Arbitrarily Anisotropic Polycrystalline Materials. Journal of Elasticity 134 (1), 1–38.

Fletcher, R., 1970. A new approach to variable metric algorithms. The Computer Journal 13, 317–322.

Fletcher, R., 2005. On the Barzilai-Borwein method. In: Qi, L., Teo, K., Yang, X. (Eds.), Optimization and Control with Applications. Springer US, Boston.

Flipon, B., Keller, C., Quey, R., Barbe, F., 2020a. A full-field crystalplasticity analysis of bimodal polycrystals. International Journal of Solids and Structures 184, 178–192.

Flipon, B., Keller, C., Quey, R., Barbe, F., 2020b. A full-field crystalplasticity analysis of bimodal polycrystals. International Journal of Solids and Structures 184, 178–192.

François, D., 1989. The Influence of the Microstructure on Fatigue. In: Moura Branco, C., Guerra Rosa, L. (Eds.), Advances in Fatigue Science and Technology. Springer Netherlands, Dordrecht, pp. 23–76.

Fraunhofer ITWM, 2021. FeelMath. Accessed: 01. July 2021.

URL https://www.itwm.fraunhofer.de/de/abteilungen/ sms/produkte-und-leistungen/feelmath.html

Frazier, P. I., Wang, J., 2016. Bayesian Optimization for Materials Design. In: Lookman, T., Alexander, F. J., Rajan, K. (Eds.), Information Science for Materials Discovery and Design. Springer, Cham, pp. 45–75.

Frederick, C. O., Armstrong, P. J., 2007. A mathematical representation of the multiaxial Bauschinger effect. Materials at High Temperatures 24 (1), 1–26.

Fritzen, F., Böhlke, T., Schnack, E., 2009. Periodic three-dimensional mesh generation for crystalline aggregates based on Voronoi tessellations. Computational Mechanics 43 (5), 701–713.

Fu, A., Zhang, J., Boyd, S., 2020. Anderson Accelerated Douglas-Rachford Splitting. SIAM Journal on Scientific Computing 42 (6), A3560– A3583.

Fullwood, D. T., Niezgoda, S. R., Adams, B. L., Kalidindi, S. R., 2010. Microstructure sensitive design for performance optimization. Progress in Materials Science 55 (6), 477–562.

Gajek, S., Schneider, M., Böhlke, T., 2021. An FE–DMN method for the multiscale analysis of short fiber reinforced plastic components. Computer Methods in Applied Mechanics and Engineering 384, 113952.

Gélébart, L., Mondon-Cancel, R., 2013. Non-linear extension of FFTbased methods accelerated by conjugate gradients to evaluate the mechanical behavior of composite materials. Computational Materials Science 77, 430 – 439.

Gel'fand, I. M., Minlos, R. A., Shapiro, Z. Y., 2018. Representations of the Rotation and Lorentz Groups and their Applications. Courier Dover Publications.

Gianola, D. S., Eberl, C., 2009. Micro-and nanoscale tensile testing of materials. JOM 61 (3), 24.

Gillner, K., Münstermann, S., 2017. Numerically predicted high cycle fatigue properties through representative volume elements of the microstructure. International Journal of Fatigue 105, 219–234.

Göküzüm, F. S., Nguyen, L. T. K., Keip, M.-A., 2019. A multiscale FE-FFT framework for electro-active materials at finite strains. Computational mechanics 64 (1), 63–84.

Goldfarb, D., 1970. A family of variable metric methods derived by variational means. Mathematics of Computation 24, 23–26.

Golovin, D., Solnik, B., Moitra, S., Kochanski, G., Karro, J., Sculley, D., 2017. Google Vizier: A Service for Black-Box Optimization. In: Proceedings of the 23rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining.

Görthofer, J., Schneider, M., Ospald, F., Hrymak, A., Böhlke, T., 2020. Computational homogenization of sheet molding compound composites based on high fidelity representative volume elements. Computational Materials Science 174, 109456.

Gottstein, G., 2013. Physical foundations of materials science. Springer, Heidelberg.

GPy, 2012. GPy: A Gaussian process framework in python. http:// github.com/SheffieldML/GPy.

Graf, M., Kuntz, M., Autenrieth, H., Diewald, F., Müller, R., 2021a. Phase field modeling and simulation of the evolution of twelve crystallographic martensite variants in austenitic parent grains. PAMM 20 (1), e202000121.

Graf, M., Kuntz, M., Autenrieth, H., Diewald, F., Müller, R., 2021b. Simulation of martensitic microstructures in a low-alloy steel. Archive of Applied Mechanics 91 (4), 1641–1668.

Griewank, A., Walther, A., 2008. Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation, 2nd Edition. SIAM, Philadelphia.

Groeber, M. A., Haley, B. K., Uchic, M. D., Dimiduk, D. M., Ghosh, S., 2006. 3D reconstruction and characterization of polycrystalline microstructures using a FIB-SEM system. Materials Characterization 57 (4), 259–273.

Groeber, M. A., Jackson, M. A., 2014. DREAM.3D: A Digital Representation Environment for the Analysis of Microstructure in 3D. Integrating Materials and Manufacturing Innovation 3 (1), 56–72.

Gross, D., Seelig, T., 2017. Fracture Mechanics: With an Introduction to Micromechanics. Springer Berlin Heidelberg, Berlin.

Guery, A., Hild, F., Latourte, F., Roux, S., 2016. Identification of crystal plasticity parameters using DIC measurements and weighted FEMU. Mechanics of Materials 100, 55–71.

Guidi, M., Adams, B. L., Onat, E. T., 1970. Tensorial Representation of the Orientation Distribution Function in Cubic Polycrystals. Textures and Microstructures 19.

Harder, J., 2001. FEM-simulation of the hardening behavior of FCC single crystals. Acta Mechanica 150 (3-4), 197–217.

Häse, F., Roch, L. M., Kreisbeck, C., Aspuru-Guzik, A., 2018. Phoenics: A Bayesian optimizer for chemistry. ACS Central Science 4 (9), 1134–1145.

Hashin, Z., Shtrikman, S., 1962a. A variational approach to the theory of the elastic behaviour of polycrystals. Journal of the Mechanics and Physics of Solids 10 (4), 343–352.

Hashin, Z., Shtrikman, S., 1962b. On some variational principles in anisotropic and nonhomogeneous elasticity. Journal of the Mechanics and Physics of Solids 10 (4), 335–342.

Hateley, J. C., Wei, H., Chen, L., 2015. Fast Methods for Computing Centroidal Voronoi Tessellations. Journal of Scientific Computing 63, 185–212.

Haupt, P., 2013. Continuum mechanics and theory of materials. Springer, New York.

Helming, K., 1996. Texturapproximation durch Modellkomponenten. Cuvillier.

Henrich, M., Pütz, F., Münstermann, S., 2020. A Novel Approach to Discrete Representative Volume Element Automation and Generation-DRAGen. Materials 13 (8), 1887.

Herrera-Solaz, V., LLorca, J., Dogan, E., Karaman, I., Segurado, J., 2014. An inverse optimization strategy to determine single crystal mechanical behavior from polycrystal tests: Application to AZ31 Mg alloy. International Journal of Plasticity 57, 1–15.

Hershey, A. V., 1954. The Elasticity of an Isotropic Aggregate of Anisotropic Cubic Crystals. Journal of Applied Mechanics 21 (3), 236– 240.

Hielscher, R., Schaeben, H., 2008. A novel pole figure inversion method: specification of the MTEX algorithm. Journal of Applied Crystallography 41 (6), 1024–1037.

Hill, R., 1952. The elastic behaviour of a crystalline aggregate. Proceedings of the Physical Society. Section A 65 (5), 349.

Hill, R., 1963. Elastic properties of reinforced solids: some theoretical principles. Journal of the Mechanics and Physics of Solids 11 (5), 357–372.

Hill, R., 1965. Continuum micro-mechanics of elastoplastic polycrystals. Journal of the Mechanics and Physics of Solids 13 (2), 89–101.

Hochhalter, J., Bomarito, G., Yeratapally, S., Leser, P., Ruggles, T., Warner, J., Leser, W., 2020. Non-deterministic Calibration of Crystal Plasticity Model Parameters. In: Ghosh, S., Woodward, C., Przybyla, C. (Eds.), Integrated Computational Materials Engineering (ICME). Springer, Cham, pp. 165–198.

Hoetzer, J., Veronika, R., Rheinheimer, W., Hoffmann, M. J., Nestler, B., 2016. Phase-field study of pore-grain boundary interaction. J Ceram Soc Jpn 124 (4), 329 – 339.

Horvath, C. D., 2021. Advanced steels for lightweight automotive structures. In: Mallick, P. K. (Ed.), Materials, design and manufacturing for lightweight vehicles. Woodhead Publishing, Sawston, pp. 39–95.

Hull, D., Bacon, D. J., 2001. Introduction to dislocations. Butterworth-Heinemann, Oxford.

Huntington, H. B., 1947. Ultrasonic measurements on single crystals. Physical Review 72 (4), 321.

Hutchinson, J. W., 1976. Bounds and self-consistent estimates for creep of polycrystalline materials. Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 348 (1652), 101–127.

Huyer, W., Neumaier, A., 1999. Global Optimization by Multilevel Coordinate Search. Journal of Global Optimization 14 (4), 331–355.

Huyer, W., Neumaier, A., 2008. SNOBFIT–Stable Noisy Optimization by Branch and Fit. ACM Transactions on Mathematical Software (TOMS) 35 (2), 1–25.

Janssens, K., 2010. An introductory review of cellular automata modeling of moving grain boundaries in polycrystalline materials. Math Comp Simul 80, 1361–1381.

Johnson, N. L., Kotz, S., Balakrishnan, N., 1994. 14: Lognormal Distributions, Continuous univariate distributions, 2nd Edition. Vol. 1 of Wiley Series in Probability and Mathematical Statistics: Applied Probability and Statistics. Wiley, New York.

Jones, D. R., Schonlau, M., Welch, W., 1998. Efficient Global Optimization of Expensive Black-Box Functions. Journal of Global Optimization 13 (4), 455–492.

Junk, M., Budday, J., Böhlke, T., 2012. On the solvability of maximum entropy moment problems in texture analysis. Mathematical Models and Methods in Applied Sciences 22 (12), 1250043.

Kabel, M., Böhlke, T., Schneider, M., 2014. Efficient fixed point and Newton-Krylov solvers for FFT-based homogenization of elasticity at large deformations. Computational Mechanics 54 (6), 1497–1514.

Kanit, T., Forest, S., Galliet, I., Mounoury, V., Jeulin, D., 2003. Determination of the size of the representative volume element for random composites: statistical and numerical approach. International Journal of Solids and Structures 40 (13-14), 3647–3679.

Kantorovich, L. V., 1948. On Newton's method for functional equations. Doklady Akademii Nauk SSSR 59, 1237–1240.

Kapoor, K., Ravi, P., Noraas, R., Park, J.-S., Venkatesh, V., Sangid, M. D., 2021. Modeling Ti–6Al–4V using crystal plasticity, calibrated with multi-scale experiments, to understand the effect of the orientation and morphology of the and phases on time dependent cyclic loading. Journal of the Mechanics and Physics of Solids 146, 104192.

Kasemer, M., Falkinger, G., Roters, F., 2020. A numerical study of the influence of crystal plasticity modeling parameters on the plastic anisotropy of rolled aluminum sheet. Modelling and Simulation in Materials Science and Engineering 28 (8), 085005.

Kawasaki, K., Nagai, T., Nakashima, K., 1989. Vertex models for twodimensional grain growth. Philos Magaz B 60 (3), 399 – 421.

Kehrer, L., Wood, J. T., Böhlke, T., 2020. Mean-field homogenization of thermoelastic material properties of a long fiber-reinforced thermoset and experimental investigation. Journal of Composite Materials 54 (25), 3777–3799.

Kitagawa, J., Mérigot, Q., Thibert, B., 2019. Convergence of a Newton algorithm for semi-discrete optimal transport. Journal of the European Mathematical Society 21 (9), 2603–2651.

Kochmann, J., Wulfinghoff, S., Reese, S., Mianroodi, J. R., Svendsen, B., 2016. Two-scale FE–FFT-and phase-field-based computational modeling of bulk microstructural evolution and macroscopic material behavior. Computer Methods in Applied Mechanics and Engineering 305, 89–110.

Kocks, U. F., Tomé, C. N., Wenk, H.-R., 2000. Texture and Anisotropy: Preferred Orientations in Polycrystals and their Effect on Materials Properties. Cambridge University Press, Cambridge.

Korte, S., Ritter, J., Jiao, C., Midgley, P. A., Clegg, W. J., 2011. Threedimensional electron backscattered diffraction analysis of deformation in MgO micropillars. Acta Materialia 59 (19), 7241–7254.

Kouznetsova, V., Brekelmans, W. A. M., Baaijens, F. P. T., 2001. An approach to micro-macro modeling of heterogeneous materials. Computational Mechanics 27, 37–48.

Krawietz, A., 1999. Parallel versus Conventional Elastoplasticity. Technische Mechanik 19 (4), 279–288.

Krill III, C., Chen, L.-Q., 2002. Computer simulation of 3-D grain growth using a phase-field model. Acta Mater 50 (12), 3059 – 3075.

Kröner, E., 1958. Berechnung der elastischen Konstanten des Vielkristalls aus den Konstanten des Einkristalls. Zeitschrift für Physik 151 (4), 504– 518.

Krupp, U., 2007. Fatigue crack propagation in metals and alloys: microstructural aspects and modelling concepts. John Wiley & Sons.

Kubis, A. J., Shiflet, G. J., Hull, R., Dunn, D. N., 2004. Focused ionbeam tomography. Metallurgical and Materials Transactions A 37 (7), 1935–1943.

Kühlbach, M., Gottstein, G., Barrales-Mora, L., 2016. A statistical ensemble cellular automaton microstructure model for primary recrystallization. Acta Mater 107, 366 – 376.

Kuhn, J., Schneider, M., Sonnweber-Ribic, P., Böhlke, T., 2020. Fast methods for computing centroidal Laguerre tessellations for prescribed volume fractions with applications to microstructure generation of polycrystalline materials. Computer Methods in Applied Mechanics and Engineering 369, 113175.

Kuhn, J., Schneider, M., Sonnweber-Ribic, P., Böhlke, T., 2022. Generating polycrystalline microstructures with prescribed tensorial texture coefficients. Computational Mechanics 70 (3), 639–659.

Kuhn, J., Spitz, J., Schneider, M., Sonnweber-Ribic, P., Böhlke, T., 2021. Identifying material parameters in crystal plasticity by Bayesian optimization. Optimization and Engineering , 1–35.

Kushner, H. J., 1964. A New Method of Locating the Maximum Point of an Arbitrary Multipeak Curve in the Presence of Noise. Journal of Basic Engineering 86, 97–106.

Lautensack, C., 2008. Fitting three-dimensional Laguerre tessellations to foam structures. Journal of Applied Statistics 35 (9), 985–995.

Lautensack, C., Zuyev, S., 2008. Random Laguerre tesselations. Advances in Applied Probability 40, 630–650.

Lebensohn, R., Rollett, A. D., 2020. Spectral methods for full-field micromechanical modelling of polycrystalline materials. Computational Materials Science 173, 109336.

Lebensohn, R. A., Tomé, C. N., neda, P. P. C., 2007. Self-consistent modelling of the mechanical behaviour of viscoplastic polycrystals incorporating intragranular field fluctuations. Philosophical Magazine 87 (28), 4287–4322.

Lehto, P., Romanoff, J., Remes, H., Sarikka, T., 2016. Characterisation of local grain size variation of welded structural steel. Welding in the World 60 (4), 673 – 688.

Lévy, B., 2015. A numerical algorithm for 2 semi-discrete optimal transport in 3D. ESAIM: Mathematical Modelling and Numerical Analysis 49 (6), 1693–1715.

Liu, W., Lian, J., Aravas, N., Münstermann, S., 2020a. A strategy for synthetic microstructure generation and crystal plasticity parameter calibration of fine-grain-structured dual-phase steel. International Journal of Plasticity 126, 102614.

Liu, W., Lian, J., Aravas, N., Münstermann, S., 2020b. A strategy for synthetic microstructure generation and crystal plasticity parameter calibration of fine-grain-structured dual-phase steel. International Journal of Plasticity 126, 102614.

Liu, W. B., Dai, Y. H., 2001. Minimization algorithms based on supervisor and searcher cooperation: I - faster and robust gradient algorithms for minimization problems with stronger noises. Journal of Optimization Theory and Applications 111, 359–379.

Liu, Y., Wang, W., Lévy, B., Sun, F., Yan, D.-M., Lu, L., Yang, C., 2016. On Centroidal Voronoi Tessellation – Energy Smoothness and Fast Computation. ACM Transaction on Graphics (TOG) 28 (4), article 244.

Lizotte, D. J., 2008. Practical Bayesian Optimization. Ph.D. thesis, University of Alberta.

Lobos, M., Böhlke, T., 2015. Materials design for the anisotropic linear elastic properties of textured cubic crystal aggregates using zeroth-, first-and second-order bounds. International Journal of Mechanics and Materials in Design 11 (1), 59–78.

Luther, T., Könke, C., 2009. Polycrystal models for the analysis of intergranular crack growth in metallic materials. Engineering Fracture Mechanics 76 (15), 2332–2343.

Lyckegaard, A., Lauridsen, E. M., Ludwig, W., Fonda, R. W., Poulsen, H. F., 2011. On the use of Laguerre tessellations for representations of 3D grain structures. Advanced Engineering Materials 13 (3), 165–170.

Mäkinen, J., 2008. Rotation manifold SO(3) and its tangential vectors. Computational Mechanics 42 (6), 907–919.

Malitsky, Y., Mishchenko, K., 2019. Adaptive gradient descent without descent. arXiv e-prints 1910.09529, 1–20.

Matérn, B., 2013. Spatial Variation. Vol. 36. Springer, New York.

Matouš, K., Geers, M. G. D., Kouznetsova, V. G., Gillman, A., 2017. A review of predictive nonlinear theories for multiscale modeling of heterogeneous materials. Journal of Computational Physics 330, 192–220.

McKay, M. D., Beckman, R. J., Conover, W. J., 2000. A Comparison of Three Methods for Selecting Values of Input Variables in the Analysis of Output From a Computer Code. Technometrics 42 (1), 55–61.

Melchior, M. A., Delannay, L., 2006. A texture discretization technique adapted to polycrystalline aggregates with non-uniform grain size. Computational Materials Science 37 (4), 557–564.

Mérigot, Q., 2011. A Multiscale Approach to Optimal Transport. Computer Graphics Forum 30, 1583–1592.

Michiuchi, M., Nambu, S., Ishimoto, Y., Inoue, J., Koseki, T., 2009. Relationship between local deformation behavior and crystallographic features of as-quenched lath martensite during uniaxial tensile deformation. Acta Materialia 57 (18), 5283–5291.

Miehe, C., Schotte, J., Schröder, J., 1999. Computational micro–macro transitions and overall moduli in the analysis of polycrystals at large strains. Computational Materials Science 16 (1-4), 372–382.

Mockus, J., 1994. Application of Bayesian approach to numerical methods of global and stochastic optimization. Journal of Global Optimization 4 (4), 347–365.

Mockus, J., Tiesis, V., Zilinskas, A., 1978. The application of Bayesian methods for seeking the extremum. In: Dixon, L. C. W., Szegö, G. P. (Eds.), Towards Global Optimisation 2. North Holland, Amsterdam, pp. 117–129.

Molinari, A., Canova, G. R., Ahzi, S., 1987. A self consistent approach of the large deformation polycrystal viscoplasticity. Acta Metallurgica 35 (12), 2983–2994.

Morawiec, A., 2003. Orientations and Rotations. Springer.

Moulinec, H., Suquet, P., 1994. A fast numerical method for computing the linear and nonlinear mechanical properties of composites. Comptes rendus de l'Académie des sciences. Série II. Mécanique, physique, chimie, astronomie. .

Moulinec, H., Suquet, P., 1998. A numerical method for computing the overall response of nonlinear composites with complex microstructure. Computer methods in applied mechanics and engineering 157 (1-2), 69–94.

MTex, 2017. MTex Documentation. https://mtextoolbox.github.io/ODFExport.html.

Müller, V., Kabel, M., Andrä, H., Böhlke, T., 2015. Homogenization of linear elastic properties of short-fiber reinforced composites–A comparison of mean field and voxel-based methods. International Journal of Solids and Structures 67, 56–70.

Nabarro, F. R. N., 1947. Dislocations in a simple cubic lattice. Proceedings of the Physical Society (1926-1948) 59 (2), 256.

Nakamura, T., Suresh, S., 1993. Effects of thermal residual stresses and fiber packing on deformation of metal-matrix composites. Acta Metallurgica et Materialia 41 (6), 1665–1681.

Natkowski, E., , Sonnweber-Ribic, P., Münstermann, S., 2021a. Determination of fatigue lifetimes with a micromechanical short crack model for the high-strength steel SAE 4150. International Journal of Fatigue , 106621.

Natkowski, E., Durmaz, A. R., Sonnweber-Ribic, P., Münstermann, S., 2021b. Fatigue lifetime prediction with a validated micromechanical short crack model for the ferritic steel EN 1.4003. International Journal of Fatigue , 106418.

Nemat-Nasser, S., Hori, M., 1993. Micromechanics: overall properties of heterogeneous materials. Elsevier, Amsterdam.

Nesterov, Y., 2004. Introductory Lectures on Convex Optimization: A Basic Course. Mathematics and its applications. Kluwer Academic Publishers, Springer US.

Nestler, B., 2005. A 3D parallel simulator for crystal growth and solidification in complex alloy systems. J Cryst Growth 275 (1), 273–278.

Nocedal, J., 1980. Updating Quasi-Newton Matrices With Limited Storage. Mathematics of Computation 35 (151), 773–782.

Nocedal, J., Wright, S. J., 1999. Numerical Optimization. Springer, New York.

Nolze, G., Hielscher, R., 2016. Orientations–perfectly colored. Journal of Applied Crystallography 49 (5), 1786–1802.

Ohno, N., Wang, J.-D., 1993. Kinematic hardening rules with critical state of dynamic recovery, part I: formulation and basic features for ratchetting behavior. International Journal of Plasticity 9 (3), 375–390.

Orowan, E., 1934. Zur kristallplastizität. i. Zeitschrift für Physik 89 (9-10), 605–613.

Pagan, D. C., Shade, P. A., Barton, N. R., Park, J.-S., Kenesei, P., Menasche, D. B., Bernier, J. V., 2017. Modeling slip system strength evolution in Ti-7Al informed by in-situ grain stress measurements. Acta Materialia 128, 406–417.

Peierls, R., 1940. The size of a dislocation. Proceedings of the Physical Society (1926-1948) 52 (1), 34.

Peirce, D., Asaro, R. J., Needleman, A., 1983. Material rate dependence and localized deformation in crystalline solids. Acta Metallurgica 31 (12), 1951–1976.

Petrich, L., Stanˇek, J., Wang, M., Westhoff, D., Heller, L., Šittner, P., Krill III, C. E., Beneš, V., Schmidt, V., 2019. Reconstruction of Grains in Polycrystalline Materials From Incomplete Data Using Laguerre Tessellations. Microscopy and Microanalysis 25, 743–752.

Picheny, V., Gramacy, R. B., Wild, S., Le Digabel, S., 2016. Bayesian optimization under mixed constraints with a slack-variable augmented Lagrangian. In: Lee, D., Sugiyama, M., Luxburg, U., Guyon, I., Garnett, R. (Eds.), Advances in Neural Information Processing Systems. Vol. 29. Curran Associates, Inc.

Polanyi, M., 1934. Über eine Art Gitterstörung, die einen Kristall plastisch machen könnte. Zeitschrift für Physik 89 (9), 660–664.

Polyak, B. T., 1987. Introduction to Optimization. Optimization Software, Inc., New York.

Prager, W., 1949. Recent Developments in the Mathematical Theory of Plasticity. Journal of Applied Physics 20 (3), 235–241.

Prasad, M. R. G., Vajragupta, N., Hartmaier, A., 2019. Kanapy: A Python package for generating complex synthetic polycrystalline microstructures. Journal of Open Source Software 4 (43), 1732.

Prithivirajan, V., Sangid, M. D., 2018. The role of defects and critical pore size analysis in the fatigue response of additively manufactured in718 via crystal plasticity. Materials & Design 150, 139–153.

Pütz, F., Henrich, M., Fehlemann, N., Roth, A., Münstermann, S., 2020. Generating Input Data for Microstructure Modelling: A Deep Learning Approach Using Generative Adversarial Networks. Materials 13 (19), 4236.

Quey, R., Dawson, P., Barbe, F., 2011. Large-scale 3D random polycrystals for the finite element method: Generation, meshing and remeshing. Computer Methods in Applied Mechanics and Engineering 200 (17-20), 1729–1745.

Quey, R., Renversade, L., 2018. Optimal polyhedral description of 3D polycrystals: Method and application to statistical and synchrotron X-ray diffraction data. Computer Methods in Applied Mechanics and Engineering 330, 308–333.

Quey, R., Villani, A., Maurice, C., 2018. Nearly uniform sampling of crystal orientations. Journal of Applied Crystallography 51 (4), 1162– 1173.

Radakrishnan, B., Zacharia, T., 1995. Monte Carlo simulation of grain boundary pinning in the weld heat-affected zone. Metall Mater Trans A 26 (8), 2123–2130.

Raponi, E., Bujny, M., Olhofer, M., Aulig, N., Boria, S., Duddeck, F., 2019. Kriging-assisted topology optimization of crash structures. Computer Methods in Applied Mechanics and Engineering 348, 730–752.

Reuss, A., 1929. Berechnung der Fließgrenze von Mischkristallen auf Grund der Plastizitätsbedingung für Einkristalle. ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik 9 (1), 49–58.

Richard, H. A., Bürgel, R., Riemer, A., 2014. Zyklische Belastung. In: Bürgel, R., Richard, H. A., Riemer, A. (Eds.), Werkstoffmechanik. Springer Vieweg, Wiesbaden, pp. 92–115.

Rieger, F., Böhlke, T., 2015. Microstructure based prediction and homogenization of the strain hardening behavior of dual-phase steel. Archive of Applied Mechanics 85 (9), 1439–1458.

Rios, L. M., Sahinidis, N. V., 2013. Derivative-free optimization: a review of algorithms and comparison of software implementations. Journal of Global Optimization 56 (3), 1247–1293.

Roe, R.-J., 1965. Description of crystallite orientation in polycrystalline materials. III. General solution to pole figure inversion. Journal of Applied Physics 36 (6), 2024–2031.

Roters, F., Eisenlohr, P., Hantcherli, L., Tjahjanto, D., Bieler, T., Raabe, D., 2010. Overview of constitutive laws, kinematics, homogenization and multiscale methods in crystal plasticity finite-element modeling: theory, experiments, applications. Acta Materialia 58, 1152–1211.

Rovinelli, A., Proudhon, H., Lebensohn, R. A., Sangid, M. D., 2020. Assessing the reliability of fast Fourier transform-based crystal plasticity simulations of a polycrystalline material near a crack tip. International Journal of Solids and Structures 184, 153–166.

Rovinelli, A., Sangid, M. D., Proudhon, H., Guilhem, Y., Lebensohn, R. A., Ludwig, W., 2018. Predicting the 3D fatigue crack growth rate of small cracks using multimodal data via Bayesian networks: In-situ experiments and crystal plasticity simulations. Journal of the Mechanics and Physics of Solids 115, 208–229.

Rycroft, C. H., 2009. Voro++: A three-dimensional Voronoi cell library in C++. Chaos 19, 041111.

Saad, Y., Schultz, M. H., 1986. GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM Journal on Scientific Computing 7, 856–869.

Saluja, R., Narayanan, R., Das, S., 2012. Cellular automata finite element (CAFE) model to predict the forming of friction stir welded blanks. Comput Mat Sci 58, 87–100.

Schäfer, B. J., Song, X., Sonnweber-Ribic, P., ul Hassan, H., Hartmaier, A., 2019a. Micromechanical Modelling of the Cyclic Deformation Behavior of Martensitic SAE - A Comparison of Different Kinematic Hardening Models. Metals 9 (3), 368.

Schäfer, B. J., Sonnweber-Ribic, P., ul Hassan, H., Hartmaier, A., 2019b. Micromechanical modeling of fatigue crack nucleation around nonmetallic inclusions in martensitic high-strength steels. Metals 9 (12), 1258.

Schemmann, M., Görthofer, J., Seelig, T., Böhlke, A. H. T., 2018. Anisotropic meanfield modeling of debonding and matrix damage in SMC composites. Composites Science and Technology 161, 143–158.

Schneider, M., 2019. On the Barzilai-Borwein basic scheme in FFT-based computational homogenization. International Journal for Numerical Methods in Engineering 118, 482–494.

Schneider, M., 2021. A review of nonlinear FFT-based computational homogenization methods. Acta Mechanica 232, 2051–2100.

Schneider, M., Josien, M., Otto, F., 2021. Representative volume elements for matrix-inclusion composites–a computational study on periodizing the ensemble. arXiv preprint arXiv:2103.07627 .

Sedighiani, K., Diehl, M., Traka, K., Roters, F., Sietsma, J., Raabe, D., 2020. An efficient and robust approach to determine material parameters of crystal plasticity constitutive laws from macro-scale stress–strain curves. International Journal of Plasticity , 102779.

Shah, A., Wilson, A., Ghahramani, Z., 2014. Student-t Processes as Alternatives to Gaussian Processes. In: Kaski, S., Corander, J. (Eds.), Proceedings of the Seventeenth International Conference on Artificial Intelligence and Statistics. Proceedings of Machine Learning Research.

Shahmardani, M., Vajragupta, N., Hartmaier, A., 2020. Robust Optimization Scheme for Inverse Method for Crystal Plasticity Model Parametrization. Materials 13 (3), 735.

Shanno, D. F., 1970. Conditioning of quasi-Newton methods for function minimization. Mathematics of Computation 24, 647–650.

Shantraj, P., Eisenlohr, P., Diehl, M., Roters, F., 2015. Numerically robust spectral methods for crystal plasticity simulations of heterogeneous materials. International Journal of Plasticity 66, 31–45.

Shenoy, M., Tjiptowidjojo, Y., McDowell, D., 2008. Microstructuresensitive modeling of polycrystalline IN 100. International Journal of Plasticity 24 (10), 1694–1730.

Simo, J. C., Hughes, T. J. R., 2006. Computational Inelasticity. Vol. 7. Springer, New York.

Smit, R. J. M., Brekelmans, W. A. M., Meijer, H. E. H., 1998. Prediction of the mechanical behavior of nonlinear heterogeneous systems by multilevel finite element modeling. Computer Methods in Applied Mechanics and Engineering 155 (1-2), 181–192.

Sobol, I. M., 1967. On the point distribution in a cube and their approximate evaluation of integrals. USSR Computational Mathematics and Mathematical Physics 7 (4), 86–112.

Spahn, J., Andrä, H., Kabel, M., Müller, R., 2014. A multiscale approach for modeling progressive damage of composite materials using fast Fourier transforms. Computer Methods in Applied Mechanics and Engineering 268, 871–883.

Spettl, A., Wertz, T., Krill III, C. E., Schmidt, V., 2014. Parametric Representation of 3D Grain Ensembles in Polycrystalline Microstructures. Journal of Statistical Physics 154, 913–928.

Spowart, J. E., Mullens, H. E., Puchalla, B. T., 2003. Collecting and analyzing microstructures in three dimensions: a fully automated approach. JOM 55 (10), 35–37.

Srinivas, N., Krause, A., Kakade, S. M., Seeger, M., 2010. Gaussian Process Optimization in the Bandit Setting: No Regret and Experimental Design. In: Fürnkranz, J., Joachims, T. (Eds.), Proceedings of the 27th International Conference on International Conference on Machine Learning. Omnipress, Madison.

Stein, M. L., 2012. Interpolation of Spatial Data: Some Theory for Kriging. Springer, New York.

Stewart, G. W., 1980. The Efficient Generation of Random Orthogonal Matrices with an Application to Condition Estimators. SIAM Journal on Numerical Analysis 17 (3), 403–409.

Student, 1908. The probable error of a mean. Biometrika , 1–25.

Tallman, A. E., Swiler, L. P., Wang, Y., McDowell, D. L., 2020. Uncertainty propagation in reduced order models based on crystal plasticity. Computer Methods in Applied Mechanics and Engineering 365, 113009.

Taylor, C. J., Kriegman, D. J., 1994. Minimization on the Lie group SO(3) and related manifolds. Yale University 16 (155), 6.

Taylor, G. I., 1934. The mechanism of plastic deformation of crystals. Part I. – Theoretical. Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 145 (855), 362–387.

Teferra, K., Rowenhorst, D. J., 2018. Direct parameter estimation for generalised balanced power diagrams. Philosophical Magazine Letters 98 (2), 79–87.

Torquato, S., 2002. Random Heterogeneous Materials: Microstructure and Macroscopic Properties. Springer, New York.

Toth, A., Kelley, C., 2015. Convergence analysis for Anderson acceleration. SIAM Journal on Numerical Analysis 53 (2), 805–819.

Tóth, L. S., Van Houtte, P., 1970. Discretization Techniques for Orientation Distribution Functions. Textures and Microstructures 19.

Tran, A., Wildey, T., 2021. Solving Stochastic Inverse Problems for Property–Structure Linkages Using Data-Consistent Inversion and Machine Learning. JOM 73 (1), 72–89.

Truesdell, C., Noll, W., 2004. The non-linear field theories of mechanics. In: Antman, S. S. (Ed.), The Non-Linear Field Theories of Mechanics. Springer Berlin, Heidelberg, Heidelberg, pp. 1–579.

Tu, X., Shahba, A., Shen, J., Ghosh, S., 2019. Microstructure and property based statistically equivalent RVEs for polycrystalline-polyphase aluminum alloys. International Journal of Plasticity 115, 268–292.

Tuegel, E. J., Ingraffea, A. R., Eason, T. G., Spottswood, S. M., 2011. Reengineering Aircraft Structural Life Prediction Using a Digital Twin. International Journal of Aerospace Engineering 2011.

Vajragupta, N., Maassen, S., Clausmeyer, T., Brands, D., Schröder, J., Hartmaier, A., 2020. Micromechanical Modeling of DP600 steel: From Microstructure to The Sheet Metal Forming Process. Procedia Manufacturing 47, 1540–1547.

Vajragupta, N., Wechsuwanmanee, P., Lian, J., Sharaf, M., Münstermann, S., Ma, A., Hartmaier, A., Bleck, W., 2014. The modeling scheme to evaluate the influence of microstructure features on microcrack formation of DP-steel: The artificial microstructure model and its application to predict the strain hardening behavior. Computational materials science 94, 198–213.

van der Walt, S., Colbert, S. C., Varoquaux, G., 2011. The NumPy array: a structure for efficient numerical computation. Computing in Science & Engineering 13 (2), 22–30.

Virtanen, P., Gommers, R., et al., T. E. O., 2020. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nature methods 17 (3), 261–272.

Voigt, W., 1889. Ueber die Beziehung zwischen den beiden Elasticitätsconstanten isotroper Körper. Annalen der physik 274 (12), 573–587.

Šilhavý, M., 1997. The Mechanics and Thermodynamics of Continuous Media. Springer, Heidelberg.

Wang, Z., Gehring, C., Kohli, P., Jegelka, S., 2018. Batched Large-scale Bayesian Optimization in High-dimensional Spaces. In: Storkey, A., Perez-Cruz, F. (Eds.), Proceedings of the Twenty-First International Conference on Artificial Intelligence and Statistics. Vol. 84. Proceedings of Machine Learning Research.

Wassermann, G., Grewen, J., 2013. Texturen metallischer Werkstoffe. Springer-Verlag.

Weygand, D., Bréchet, Y., Lépinoux, J., Gust, W., 1999. Three-dimensional grain growth: a vertex dynamics simulation. Philos Magaz B 79 (5), 703–716.

Wicht, D., Schneider, M., Böhlke, T., 2020a. An efficient solution scheme for small-strain crystal-elasto-viscoplasticity in a dual framework. Computer Methods in Applied Mechanics and Engineering 358, 112611.

Wicht, D., Schneider, M., Böhlke, T., 2020b. On Quasi-Newton methods in FFT-based micromechanics. International Journal for Numerical Methods in Engineering online, 1–34.

Williams, C. K. I., Rasmussen, C. E., 2006. Gaussian Processes for Machine Learning. Vol. 2. The MIT press, Cambridge.

Wilson, J. T., Hutter, F., Deisenroth, M. P., 2018. Maximizing acquisition functions for Bayesian optimization. In: 32nd Conference on Neural Information Processing Systems (NeurIPS 2018).

Wu, B., Vajragupta, N., Lian, J., Hangen, U., Wechsuwanmanee, P., Münstermann, S., 2017. Prediction of plasticity and damage initiation behaviour of C45E+ N steel by micromechanical modelling. Materials & Design 121, 154–166.

Wulfinghoff, S., Bayerschen, E., Böhlke, T., 2013. A gradient plasticity grain boundary yield theory. International Journal of Plasticity 51, 33–46.

Yang, S., Dirrenberger, J., Monteiro, E., Ranc, N., 2019. Representative volume element size determination for viscoplastic properties in polycrystalline materials. International Journal of Solids and Structures 158, 210–219.

Zaefferer, S., Wright, S. I., Raabe, D., 2008. Three-dimensional orientation microscopy in a focused ion beam-scanning electron microscope: a new dimension of microstructure characterization. Metallurgical and Materials Transactions A 39 (2), 374–389.

Zeman, J., Vondˇrejc, J., Novák, J., Marek, I., 2010. Accelerating a FFT-based solver for numerical homogenization of periodic media by conjugate gradients. Journal of Computational Physics 229 (21), 8065–8071.

Zheng, Q.-S., Zou, Y.-B., 2001a. Orientation Distribution Functions for Microstructures of Heterogeneous Materials (I)–Directional Distribution Functions and Irreducible Tensors. Applied Mathematics and Mechanics 22 (8), 865–884.

Zheng, Q.-S., Zou, Y.-B., 2001b. Orientation Distribution Functions for Microstructures of Heterogeneous Materials (II)–Crystal Distribution Functions and Irreducible Tensors Restricted by Various Material Symmetries. Applied Mathematics and Mechanics 22 (8), 885–903.

Zhuang, X., Wang, Q., Zhu, H., 2015. A 3D computational homogenization model for porous material and parameters identification. Computational Materials Science 96, 536–548.

#### **Schriftenreihe Kontinuumsmechanik im Maschinenbau Karlsruher Institut für Technologie (KIT) (ISSN 2192-693X)**


Die Bände sind unter www.ksp.kit.edu als PDF frei verfügbar oder als Druckausgabe bestellbar.


Die Bände sind unter www.ksp.kit.edu als PDF frei verfügbar oder als Druckausgabe bestellbar.


#### Juliane Lang **Thermomechanical Modeling and Experimental Characterization of Sheet Molding Compound Composites.**  ISBN 978-3-7315-1232-5 **Band 22**

Die Bände sind unter www.ksp.kit.edu als PDF frei verfügbar oder als Druckausgabe bestellbar.


To capture the influence of the microstructure on the cyclic mechanical behavior of polycrystalline metals in a resource-efficient way, computational homogenization methods solve equations on computational cells. These reflect the essential characteristics of the microstructure. In this work we propose to use modern solvers to compute representations of polycrystalline materials in terms of Laguerre tessellations with prescribed volume fractions. We assess the performance of these solvers in terms of the number of iterations and overall run-time needed to compute the Laguerre cells. To assign a crystallographic orientation to each cell, we propose an optimization scheme based on tensorial texture coefficients, a measure used to quantify the anisotropy of polycrystals. We investigate the capability of the method to reproduce isotropic and anisotropic mechanical behavior with microstructures featuring different grain size distributions. Describing the deformation behavior on a microscopic scale requires identifying suitable parameters for the material model at hand, i.e., single crystal plasticity. As the ground truth we use macroscopic experiments and compare them with computational results, leading to an inverse optimization problem. We propose to use a Bayesian optimization strategy, investigate its performance as well as its robustness and compare it to other algorithms.

ISSN 2192-693X ISBN 978-3-7315-1272-1

Gedruckt auf FSC-zertifiziertem Papier **25**

J. Kuhn **Micromechanical modeling of polycrystalline metals**